tcb@breakable [1]∫⋯∫

#  Efficient Bilevel Source Mask Optimization

Guojin Chen Chinese University of Hong Kong<EMAIL_ADDRESS>,
Hongquan He Shanghai Tech. University<EMAIL_ADDRESS>,  Peng Xu Chinese
University of Hong Kong<EMAIL_ADDRESS>,  Hao Geng Shanghai Tech.
University<EMAIL_ADDRESS>and  Bei Yu Chinese University of Hong
Kong<EMAIL_ADDRESS>
###### Abstract.

Resolution Enhancement Techniques (RETs) are critical to meet the demands of
advanced technology nodes. Among RETs, Source Mask Optimization (SMO) is
pivotal, concurrently optimizing both the source and the mask to expand the
process window. Traditional SMO methods, however, are limited by sequential
and alternating optimizations, leading to extended runtimes without
performance guarantees. This paper introduces a unified SMO framework
utilizing the accelerated Abbe forward imaging to enhance precision and
efficiency. Further, we propose the innovative BiSMO framework, which
reformulates SMO through a bilevel optimization approach, and present three
gradient-based methods to tackle the challenges of bilevel SMO. Our
experimental results demonstrate that BiSMO achieves a remarkable 40%
reduction in error metrics and 8$\times$ increase in runtime efficiency,
signifying a major leap forward in SMO.

††copyright: none

##  1\. Introduction

Lithography, vital for semiconductor manufacturing, advances integrated
circuit (IC) development. The semiconductor industry’s drive for
miniaturization and efficiency challenges traditional lithography in creating
finer patterns, emphasizing the need for resolution enhancement techniques
(RETs) to meet advanced semiconductor requirements. Among various RETs such as
sub-resolution assist features (SRAFs) (OPC-USP2009-Wallace, ), optical
proximity correction (OPC) (OPC-DAC2014-Gao, ), and source mask optimization
(SMO), SMO stands out due to its broader solution space. SMO uniquely
optimizes both illumination source and mask pattern, ensuring pivotal
lithographic fidelity vital for advancing next-generation IC manufacturing.

As shown in Figure 1(a), SMO integrates source optimization (SO), mask
optimization (MO), and their iterative optimization refinement. The efficiency
and performance of SMO are primarily influenced by two factors: the imaging
model and the optimization strategy. Central to both SO and MO are the forward
imaging models: Abbe’s (Fuehner2009abbe, ) and Hopkins’ (cobb1995sum, ), each
with unique computational characteristics. Abbe’s model, celebrated for its
precision, demands extensive computation through the summation of intensities
from discrete source points. Conversely, Hopkins’ model, employing truncated
singular value decomposition (SVD), reduces computational load but is
unsuitable for SO due to its loss of source information. Beyond isolated SO
and MO, the strategy employed in their combined optimization significantly
influences the optimization trajectory and solution space, thereby affecting
SMO outcomes.

Within MO, Hopkins’ model is foundational to various methods. This includes
MOSAIC (OPC-DAC2014-Gao, ) which blends design target and process window
considerations, and techniques such as GAN-OPC (OPC-DAC2018-Yang, ), DAMO
(OPC-ICCAD2020-DAMO, ), Neural-ILT (OPC-ICCAD2020-NeuralILT, ), and DevelSet
(OPC-ICCAD2021-DevelSet, ) that employ deep neural networks for enhanced
performance. Hardware acceleration strategies like those in GPU-LS (OPC-
DATE21-Yu, ) and Multi-ILT (multiilt2023sun, ) leverage GPU and parallel
acceleration to expedite Hopkins-based MO tasks. To the best of our knowledge,
the potential of Abbe-based MO optimization and acceleration remains
unexplored.

For SO, the impracticality of Hopkins’ model necessitates exclusive reliance
on Abbe’s. Previous SO strategies have employed compressive sensing
(wang2020csfasticc, ) and sampling-based methods (sampleICC2022sun, ) for
complexity reduction. However, the benefit is limited. By contrast, MO using
the Hopkins model with GPU or DNN acceleration achieve average optimization
times of five seconds, while SO dependent on Abbe’s model typically exceed 30
minutes. Due to its computational intensity, accelerating Abbe’s model is
essential.

Figure 1. (a) The forward lithography and SMO process. (b) Bilevel SMO with
upper-level MO and lower-level SO.

In SMO, alongside the imaging model, the co-optimization strategy is a
critical determinant of performance. Since the inception of SMO in early
2010s, the alternating minimization (AM) strategy has been a dominant approach
in SMO (wang2020csfasticc, ; ding2020gradient, ; sampleICC2022sun, ). As shown
in Figure 2(a), AM method involves isolated iterative minimization of source
parameters across multiple SO epochs while maintaining constant mask
variables, followed by mask optimization in MO epochs with the source fixed,
repeating until convergence. However, AM strategy’s simplicity does not
guarantee effectiveness. Due to the concurrent impact of source and mask
parameters on the aerial image, AM’s localized focus can result in SMO being
confined to local minima, and the absence of global, gradient-based guidance
can prolong convergence, necessitating more iterations. These limitations
highlight the need for more advanced co-optimization strategies in SMO.

To address SMO challenges, we utilize the Abbe model for forward imaging,
leveraging its capability for concurrent SO and MO gradient computation, and
its superior precision due to avoidance of approximate decomposition. We
counterbalance the Abbe model’s intensive computational demand with GPU
acceleration. This enables us to develop an Abbe-based unified SMO framework
incorporating process window considerations, consequently enhancing SMO
outcomes, improving efficiency, and reducing process variability. Furthermore,
we innovatively reconceptualize the SMO problem as a bilevel optimization
(BLO) challenge to gain a better co-optimization strategy, as depicted in
Figure 1(b). BLO, a hierarchical mathematical program, is defined as an
optimization problem where the feasible region is constrained by another
nested optimization problem. It is widely applied in areas such as
hyperparameter optimization (lorraine2020optimizing, ), neural architecture
search (zhang2021idarts, ), and multitask or meta-learning (choe2023making, ).
We then propose three novel gradient-based bilevel SMO solutions, featuring
global perspectives achieved through MO and SO gradient fusion. These
solutions demonstrate substantial improvements in efficiency and accuracy over
the traditional SMO approaches. Our primary contributions are as follows:

  * •

We establish the first unified Abbe-based SMO framework incorporating process
window considerations, significantly accelerating Abbe imaging through
parallel computation, achieving speeds comparable to Hopkins’ method.

  * •

We pioneer the modeling of SMO as a unified bilevel framework, developing
three efficient gradient-based methods with global perspectives and improved
exploration of solution space, superseding the conventional SMO.

  * •

Our experimental results indicate that, compared to state-of-the-art (SOTA)
SMO methods (sampleICC2022sun, ), our approach reduces error metrics by
approximately 40% and achieves an eightfold increase in throughput. In
comparison to SOTA MO methods (multiilt2023sun, ), our error metrics are half
as large.

##  2\. Preliminaries

###  2.1. Lithography Simulation

In optical lithography systems, the intensity of aerial image $I(x,y)$ on the
wafer plane can be formulated via lithography theory (cobb1995sum, ) as:

(1)

$\begin{aligned}
{I}(x,y)=&\\!\multiint{6}_{-\infty}^{\infty}{J}(f,g)O(f^{\prime},g^{\prime})O^{*}(f^{\prime\prime},g^{\prime\prime}){H}(f\\!+\\!f^{\prime},g\\!+\\!g^{\prime})\\\
&{H}^{*}(f\\!+\\!f^{\prime\prime},\\!g+g^{\prime\prime})\exp(-i2\pi((f^{\prime}\\!-\\!f^{\prime\prime})x+(g^{\prime}\\!-\\!g^{\prime\prime})y))\\\
&\mathrm{d}f~{}\mathrm{d}g~{}\mathrm{d}f^{\prime}~{}\mathrm{d}g^{\prime}~{}\mathrm{d}f^{\prime\prime}~{}\mathrm{d}g^{\prime\prime},\end{aligned}$

where $\bm{J}$ is the illumination source. $\bm{H}$ is projection system’s
optical transfer function. $O(f^{\prime},g^{\prime})$ captures the binary mask
pattern $M(x,y)$’s frequency spectrum, derived via a $2D$ fast Fourier
transform (FFT), $\mathcal{F}(\cdot)$. ∗ signifies the Hermitian transpose.
$(x,y)$ denotes spatial coordinates, while $(f,g)$, $(f^{\prime},g^{\prime})$,
and $(f^{\prime\prime},g^{\prime\prime})$ refer to the frequency coordinates
of the source, mask spectrum, and its conjugate, respectively. The formation
of aerial image in Equation 1 is computed using two distinct methods: Abbe’s
and Hopkins’ method.

Abbe’s Approach: Abbe’s approach, also known as the source points integration
approach, discretizes the source space and independently computes the
contribution of each source point, subsequently summing these contributions to
form the aerial image. Regardless of the discretization technique, the source
can hence be represented as a set of source points
$\\{(f_{\sigma},g_{\sigma};j_{\sigma})\\}$, where each source point defines a
pair of spatial frequencies and its discrete magnitude $j_{\sigma}\in[0,1]$.
By setting
$\bm{A}_{(f,g)}(f^{\prime},g^{\prime})={H}(f+f^{\prime},g+g^{\prime})O(f^{\prime},g^{\prime}),$
and applying the Inverse Fast Fourier Transform (IFFT), the total intensity in
Equation 1 can be formulated in Abbe’s approach:

(2)
$\displaystyle{I}(x,y)=\sum_{\sigma}j_{\sigma}|\bm{A}_{(f_{\sigma},g_{\sigma})}(x,y)|^{2}.$

Hopkins’ Approach: Hopkins’ approach separates the calculation of source and
projection system from the processing of the mask for Equation 1. It
formulates the source and projector into the transmission cross-coefficients
($TCC$), as defined by:

(3)

$\begin{aligned}
{TCC}=\\!\iint_{-\infty}^{\infty}{J}(f,g){H}(f\\!+\\!f^{\prime}\\!,g\\!+\\!g^{\prime}){H}^{*}(f\\!+\\!f^{\prime\prime},\\!g\\!+\\!g^{\prime\prime})\mathrm{d}f\mathrm{d}g.\end{aligned}$

The Sum of Coherent Systems (SOCS) (cobb1995sum, ) provides an approximation
to the Hopkins imaging equations, simplifying the ${TCC}$ spectrum using SVD.
Due to the rapid decay of eigenvalue $\kappa_{q}$ with $q$, only the top
truncated $Q$ eigenvalues are retained. By applying the IFFT, SOCS can be
expressed in spatial domain as:

(4) ${I}(x,y)=\sum_{q=1}^{Q}\kappa_{q}|{\phi}_{q}(x,y)\otimes{M}(x,y)|^{2},$

where ${\phi}_{q}(x,y)$ is the spatial distribution of eigenvector
$\bm{\Phi_{q}}$. Here, $\otimes$ denotes convolution and $|\cdot|$ is the
absolute operator. Hopkins’ Approach vs. Abbe’s Approach: In Equation 4,
Hopkins’ method reduces computational demands from Abbe’s
$\mathcal{O}(N_{j}^{2}\cdot N_{m}^{4})$ to $\mathcal{O}(Q\cdot N_{m}^{4})$,
with $Q<N_{j}^{2}$ for source $\bm{J}\in\mathbb{R}^{N_{j}\times N_{j}}$ and
mask $\bm{M}\in\mathbb{R}^{N_{m}\times N_{m}}$. Under identical optical
conditions, Hopkins’ method outperforms Abbe’s in speed, leading to its
preference in various MO algorithms (OPC-DAC2014-Gao, ; OPC-DAC2018-Yang, ;
OPC-ICCAD2020-DAMO, ; OPC-ICCAD2020-NeuralILT, ; OPC-ICCAD2021-DevelSet, ;
OPC-DATE21-Yu, ; multiilt2023sun, ). However, Hopkins’ reliance on truncated
SVD, as shown in Equation 4, prevents SO due to the inability to calculate
source gradients. In contrast, Abbe’s method (Equation 2) inherently suits SO
by summing the impacts of all source points to form the aerial image.
Moreover, Abbe’s richer source information enhances lithography precision,
thus improving MO outcomes. Consequently, Abbe’s method is essential for SO,
higher MO precision, and indispensable for gradient-based bilevel SMO.

###  2.2. Evaluation Metrics

######  Definition 1 (Squared $L_{2}$ Error (L2)).

Given target pattern $\bm{Z}_{t}$ and resist image under nominal process
condition $\bm{Z}$, the squared $L_{2}$ error is calculated as
$\|\bm{Z}-\bm{Z}_{t}\|_{2}^{2}$.

######  Definition 2 (Process Variation Band (PVB)).

PVB (OPC-ICCAD2013-Banerjee, ) is used in manufacturing to represent the
expected range of variation in a production process. PVB denotes the
$\operatorname{XOR}$ area between the resist images $\bm{Z}_{\min}$ and
$\bm{Z}_{\max}$ under the min and max process conditions.

######  Definition 3 (Edge Placement Error (EPE)).

EPE (OPC-ICCAD2013-Banerjee, ) refers to the deviation between the intended
position of a feature on a wafer and its actual position after lithography.

##  3\. Algorithm

In this section, we detail the core contributions of our study: the
development and acceleration of an Abbe-based SMO framework with process
window considerations in Section 3.1. Furthermore, in Section 3.2, we offer a
thorough derivation and comparison of gradient-based bilevel SMO, where the
upper-level MO gains a global perspective by gradient fusion, enabling more
informed decision-makeing by solving the lower-level SO.

###  3.1. Abbe-based Unified SMO Framework

The Hopkins model, hindered by the need for frequent, inefficient, and non-
differentiable SVD truncation of the TCC, limits effective gradient-based co-
optimization of source and mask. Addressing this, we introduce a unified SMO
framework employing the Abbe model, which facilitates efficient, joint SMO
without the TCC’s burdensome processing. The framework is further enhanced by
incorporating parallel computing techniques for acceleration.

Abbe-based unified SMO: Utilizing freeform illumination, the pixelated source
point is denoted as $J(f,g)\in[0,1]$, while binary mask values $M(x,y)$ are
either 0 or 1. To render the SMO framework differentiable, we introduce
optimization parameters $\bm{\theta}_{J}$ for Source $\bm{J}$ and
$\bm{\theta}_{M}$ for Mask $\bm{M}$, where
$\bm{\theta}_{J}\in\mathbb{R}^{N_{j}\times N_{j}}$,
$\bm{\theta}_{M}\in\mathbb{R}^{N_{m}\times N_{m}}$, and both parameters can
assume any real value. Here, $N_{j}$ and $N_{m}$ represent the dimensions of
the source and mask. Appropriate activation and initialization enable deriving
source and mask from these parameters. The Sigmoid function
$\sigma(x)=1/\left(1+\exp(-x)\right)$ is employed for both grayscale source
and binary mask, as listed in Table 1.

Table 1. The activation and initialization for Abbe-imaging.

| Activation | Initialization
---|---|---
Mask $\bm{M}$ | $\bm{M}=\sigma(\alpha_{m}\cdot\bm{\theta}_{M})$ | $\bm{\theta}_{M}(x,y)={m_{0}},\text{ if }M_{0}(x,y)=1;\text{else }-\\!m_{0}.$
Source $\bm{J}$ | $\bm{J}=\sigma(\alpha_{j}\cdot\bm{\theta}_{J})$ | $\bm{\theta}_{J}(f,g)={j_{0}},\text{ if }J_{0}(f,g)=1;\text{ else }-\\!j_{0}.$

In Table 1, $\alpha_{m}$ and $\alpha_{j}$ are the sigmoid steepness. Initial
values for $\bm{\theta}_{M}$ are assigned as either $m_{0}$ or $-m_{0}$ based
on initial mask pattern $\bm{M}_{0}$. Typically, the initial mask pattern is
the same as the binary target pattern $\bm{Z}_{t}$. This mask initialization
also facilitates SRAF generation during MO. Grayscale source $\bm{J}$ requires
careful selection of steepness $\alpha_{j}$ and initialization $j_{0}$ to
maintain its grayscale property. Hyperparameters are detailed in Section 4.
The shape of initial source pattern $\bm{J}_{0}$ is derived from parametric
templates like annular, quasar, or dipole, characterized by outer and inner
radii $\sigma_{o}$ and $\sigma_{i}$. Although the Cosine function is an
alternative source activation, its use may lead to training instability due to
gradient issues, leading us to prefer the Sigmoid function. The transfer
function $H(f,g)$ can be accurately characterized by a low-pass filter,
expressed as:

(5) $H(f,g)=\begin{cases}1,&\text{if
}\sqrt{f^{2}+g^{2}}\leq\frac{{N\\!A}}{\lambda},\\\
0,&\text{otherwise},\end{cases}$

where the cut-off frequency is determined by the projection system’s numerical
aperture ${N\\!A}$ and the illumination wavelength $\lambda$. By integrating
Equation 2, Equation 5 and Table 1, we formulate the Abbe forward imaging
$f_{abbe}$, determining the aerial image $\bm{I}$ as a function of the
optimization parameters for source $\bm{\theta}_{J}$ and mask
$\bm{\theta}_{M}$: $\bm{I}=f_{abbe}(\bm{\theta}_{J},\bm{\theta}_{M}).$ Then we
can utilize a straightforward threshold model for resist modeling. The Sigmoid
activation is also adopted to ensure a smooth transition and maintain
differentiability:

(6) $\bm{Z}=\sigma(\beta\cdot(\bm{I}-{I}_{tr})),$

where $\bm{Z}$ represents the resist pattern, ${I}_{tr}$ denotes the intensity
threshold, and $\beta$ is the steepness.

We have now established a complete Abbe forward imaging model that maps the
parameters $\bm{\theta}_{J}$ and $\bm{\theta}_{M}$ to the aerial image
$\bm{I}$ and the resist image $\bm{Z}$. To realize SMO, it is essential to
define the corresponding objective function and optimization method. We employ
the mean squared loss to quantify the discrepancy between the resist pattern
$\bm{Z}$ and the target pattern $\bm{Z}_{t}$:

(7) $\mathcal{L}_{2}=\|\bm{Z}-\bm{Z}_{t}\|^{2}.$

In alignment with (OPC-DAC2014-Gao, ) and considering a $\pm 2\%$ dose range,
we pioneer the integration of process window considerations into Abbe-based
SMO to mitigate process variation via PVB loss. By substituting
$\bm{M}_{\min}\\!=\\!d_{\min}\cdot\sigma(\alpha_{m}\cdot\bm{\theta}_{M})$ and
$\bm{M}_{\max}\\!=\\!d_{\max}\cdot\sigma(\alpha_{m}\cdot\bm{\theta}_{M})$ into
$f_{abbe}$, we obtain the resist patterns $\bm{Z}_{\min}$, $\bm{Z}_{\max}$
under minimum $d_{\min}$ and maximum $d_{\max}$ process conditions. The PVB
loss is formulated as:

(8)
$\mathcal{L}_{pvb}=\|\bm{Z}_{\max}-\bm{Z}_{t}\|^{2}+\|\bm{Z}_{\min}-\bm{Z}_{t}\|^{2}.$

Consequently, the comprehensive SMO loss $\mathcal{L}_{smo}$ is formulated as:

(9)
$\mathcal{L}_{smo}\coloneqq\mathcal{L}_{so}\coloneqq\mathcal{L}_{mo}=\gamma\mathcal{L}_{2}+\eta\mathcal{L}_{pvb},$

where $\gamma$ and $\eta$ are weighting factors for the respective loss
components. SO loss $\mathcal{L}_{so}$ and MO loss $\mathcal{L}_{mo}$ can
utilize same objective functions. The SMO problem is thus defined as:

(10)
$(\hat{\bm{\theta}_{J}},\hat{\bm{\theta}_{M}})=\underset{(\bm{\theta}_{J},\bm{\theta}_{M})}{\operatorname{argmin}}~{}\mathcal{L}_{smo}(\bm{\theta}_{J},\bm{\theta}_{M}),$

where $\hat{\bm{\theta}_{J}},\hat{\bm{\theta}_{M}}$ represent the optimal
parameter values for the source and mask, respectively.

Abbe acceleration. The primary computational bottleneck in SMO is the forward
imaging model and its gradient calculations. As noted in Equations 2 and 4,
contributions from source points can be independently calculated, making the
complexity ratio between Abbe’s and Hopkins’s models $\sigma/Q$, where
$\sigma\in[0,N_{j}^{2}]$ represents the number of effective source points
(i.e., where $j_{\sigma}>0$). Both models can be accelerated with parallel
computing, using multicore CPUs or GPUs. Theoretically, the parallel
computation time ratio is
$\lceil\frac{\sigma}{P}\rceil/\lceil\frac{Q}{P}\rceil$, where $P$ denotes the
maximum number of parallel threads and $\lceil\cdot\rceil$ is the ceiling
operator. This suggests that Abbe’s runtime can match Hopkins’ if
$P\geq\sigma$. In our implementation, GPUs are utilized for parallel
computation of each effective source point’s contribution to the aerial image,
owing to their greater thread parallelism, larger memory bandwidth, and faster
FFT or IFFT operations compared to multicore CPUs. Experimental results in
Figure 3 demonstrate that our Abbe-imaging model achieves runtime performance
comparable to Hopkins’ model.

###  3.2. Efficient Bilevel SMO

Previous alternating minimization-based SMO: Since the introduction of SMO
technology, the significant computational demands have compelled previous
methods to compromise on a simple alternating minimization (AM) strategy. As
illustrated in Algorithm 1, AM-based SMO (AM-SMO) alternates between two
minimization cycles, updating parameters $\bm{\theta}_{J}$ and
$\bm{\theta}_{M}$ sequentially. This process iterates until reaching the
specified convergence criteria for SO and MO, as shown in Figure 2(a).
However, AM-SMO has several notable drawbacks: 1) AM-SMO tends to converge to
local minima due to its narrow focus on localized aspects of SO or MO,
ignoring the global structure of the problem. 2) The convergence is often slow
because the source and mask are highly interdependent. Adjusting
$(\bm{\theta}_{M})_{k}$ as per 5 makes $(\bm{\theta}_{J})_{k}$ suboptimal,
requiring numerous iterations for stabilization. 3) The absence of global
gradient guidance complicates establishing effective early stopping criteria,
often resulting in either prolonged optimization or suboptimal convergence.

Algorithm 1 Alternating Minimization-based SMO (AM-SMO) (wang2020csfasticc, ;
sampleICC2022sun, )

1: for $k=1,2,3,\ldots$; do $\rhd$ Alternating SO & MO. for $k=1,2,3,\ldots$;
do $\rhd$ Alternating SO & MO.

2: for $k=1,2,3,\ldots$; do $\rhd$ Alternating SO & MO. for $k=1,2,3,\ldots$;
do $\rhd$ Alternating SO & MO. for $k=1,2,3,\ldots$; do $\rhd$ Alternating SO
& MO. while not converged do $\rhd$ SO iterations. for $k=1,2,3,\ldots$; do
$\rhd$ Alternating SO & MO. for $k=1,2,3,\ldots$; do $\rhd$ Alternating SO &
MO. while not converged do $\rhd$ SO iterations.

3: for $k=1,2,3,\ldots$; do $\rhd$ Alternating SO & MO. for $k=1,2,3,\ldots$;
do $\rhd$ Alternating SO & MO. while not converged do $\rhd$ SO iterations.
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed.
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed.
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed.

4:
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed.
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed.
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed. while not converged do $\rhd$ MO
iterations.
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed.
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed. while not converged do $\rhd$ MO
iterations.

5:
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed.
$(\bm{\theta}_{J})_{k}\leftarrow{\operatorname{argmin}}_{\bm{\theta}_{J}}\mathcal{L}_{so}\left(\bm{\theta}_{J},(\bm{\theta}_{M})_{k-1}\right)$;
$\rhd$ $\bm{\theta}_{M}$ is fixed. while not converged do $\rhd$ MO
iterations.
$(\bm{\theta}_{M})_{k}\leftarrow\operatorname{argmin}_{\bm{\theta}_{M}}\mathcal{L}_{mo}\left((\bm{\theta}_{J})_{k},\bm{\theta}_{M}\right)$;
$\rhd$ $\bm{\theta}_{J}$ is fixed.
$(\bm{\theta}_{M})_{k}\leftarrow\operatorname{argmin}_{\bm{\theta}_{M}}\mathcal{L}_{mo}\left((\bm{\theta}_{J})_{k},\bm{\theta}_{M}\right)$;
$\rhd$ $\bm{\theta}_{J}$ is fixed.
$(\bm{\theta}_{M})_{k}\leftarrow\operatorname{argmin}_{\bm{\theta}_{M}}\mathcal{L}_{mo}\left((\bm{\theta}_{J})_{k},\bm{\theta}_{M}\right)$;
$\rhd$ $\bm{\theta}_{J}$ is fixed. return
$(\bm{\theta}_{J})_{k},(\bm{\theta}_{M})_{k}$.

Proposed bilevel SMO: Equation 10 frames SMO as a typical multivariate
optimization problem. Yet, the AM-SMO, as detailed in Algorithm 1, often leads
to suboptimal results due to its localized focus. A gradient-based approach
with a global perspective is essential to overcome these limitations. In SMO,
SO and MO have a nested relationship. The goal is to rapidly and effectively
determine the optimal source response for each altered mask, offering MO
global gradient direction to boost overall SMO efficiency. This scenario is
particularly well-suited for a bilevel optimization approach. From a bilevel
viewpoint, the upper-level MO, constrained by the optimal solutions from the
lower-level SO, forms a dependent hierarchy, as illustrated in Figure 1(b) and
Figure 2(b). This structure allows the MO to offer a global perspective by
solving SO, guiding the optimization beyond local minima. Consequently,
Equation 10 can be reformulated in a bilevel context:

(11) $\displaystyle\begin{aligned} \min_{\bm{\theta}_{M}}&\
\mathcal{L}_{mo}(\bm{\theta}_{J}^{*}(\bm{\theta}_{M}),\bm{\theta}_{M}),&\leavevmode\hbox{}$$\halign{#\cr\vskip
0.0pt\vskip-99.0pt\cr}$$\vrule width=0.0pt,depth=0.0pt\nobreak\hskip
0.0pt\ignorespaces\leavevmode\hbox{}$$\halign{#\cr\vskip
0.0pt\vskip-99.0pt\cr}$$\vrule width=0.0pt,depth=0.0pt\nobreak\hskip
363.48921pt\ignorespaces\textit{{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>Level: MO}}\ignorespaces\\\ \textup{s.t.}&\
\bm{\theta}_{J}^{*}(\bm{\theta}_{M})={\operatorname{argmin}}_{\bm{\theta}_{J}}~{}\mathcal{L}_{so}(\bm{\theta}_{J},\bm{\theta}_{M}).&\leavevmode\hbox{}$$\halign{#\cr\vskip
0.0pt\vskip-99.0pt\cr}$$\vrule width=0.0pt,depth=0.0pt\nobreak\hskip
0.0pt\ignorespaces\leavevmode\hbox{}$$\halign{#\cr\vskip
0.0pt\vskip-99.0pt\cr}$$\vrule width=0.0pt,depth=0.0pt\nobreak\hskip
367.1559pt\ignorespaces\textit{{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>Level: SO}}\ignorespaces\end{aligned}$

In BLO, the inner and outer loops are termed lower-level and upper-level
subproblems. Here the inner loop seeks optimal source parameters
$\bm{\theta}_{J}^{*}$ for the current $\bm{\theta}_{M}$ while the outer loop
endeavors to optimize the mask parameters $\bm{\theta}_{M}$ with the best-
response source $\bm{\theta}_{J}^{*}(\bm{\theta}_{M})$. The gradient of the
outer loop in bilevel SMO, also referred to as the hypergradient, which is
derived from the fusion of gradients from the upper and lower levels, is then
calculated as:

(12)
$\nabla_{\bm{\theta}_{M}}\mathcal{L}_{mo}=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{M}}+\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\frac{\partial\bm{\theta}_{J}^{*}(\bm{\theta}_{M})}{\partial\bm{\theta}_{M}}.$

Within the Abbe-based SMO framework introduced in Section 3.1, the direct
gradient $\frac{\partial\mathcal{L}{mo}}{\partial\bm{\theta}_{M}}$ and
$\frac{\partial\mathcal{L}{mo}}{\partial\bm{\theta}_{J}}$ can be efficiently
computed. However, two principal challenges arise: (1) the precise
approximation of the SO optimal solution
$\bm{\theta}_{J}^{*}(\bm{\theta}_{M})$, and (2) differentiating the best-
response Jacobian:
$\frac{\partial\bm{\theta}_{J}^{*}(\bm{\theta}_{M})}{\partial\bm{\theta}_{M}}$.
To address the former, we approximate $\bm{\theta}_{J}^{*}$ by unrolling a few
SO gradient steps $T$, significantly reducing the computational cost.
Fortunately, extensive research (lorraine2020optimizing, ; zhang2021idarts, ;
choe2023making, ) indicates that BLO, with weight sharing, can effectively
adapt $\bm{\theta}_{J}$ to $\bm{\theta}_{J}^{*}$ with a small unrolling step
$T$. For the latter issue, we propose three methods to compute the best-
response Jacobian: bilevel SMO using finite difference (BiSMO-FD), BiSMO-NMN
utilizing Neumann series, and BiSMO-CG using conjugate gradients.

Figure 2. (a)Previous AM-SMO flow. (b)Our BiSMO flow.

####  3.2.1. BiSMO-FD:

The finite difference (FD) strategy uses a single inner SO step,
$\bm{\theta}_{J}^{*}(\bm{\theta}_{M})=\bm{\theta}_{J}-\xi\nabla_{\bm{\theta}_{J}}\mathcal{L}_{so}$,
and obtaining
$\frac{\partial\bm{\theta}_{J}^{*}(\bm{\theta}_{M})}{\partial\bm{\theta}_{M}}=-\xi\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}$,
the BiSMO-FD calculates the hypergradient as:

(13)
$\nabla_{\bm{\theta}_{M}}\mathcal{L}_{mo}^{\texttt{FD}}=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{M}}-\xi\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}},$

where $\xi$ is the inner-loop learning rate for $\bm{\theta}_{J}$. While
increasing the number of inner SO steps can lead to a more precise
approximation for $\bm{\theta}_{J}^{*}$, akin to AM-SMO, it results in a
linear increase in memory and computational load, becoming impractical due to
the need to store optimization paths and all intermediate gradients for
differentiation following the chain rule.

In contrast, using the implicit function theorem (IFT), the hypergradient can
be computed without retaining intermediate gradients, thus independent of the
optimization path and significantly reducing memory usage. The IFT-based
hypergradient for SMO can be formulated as the following lemma.

######  Lemma 1.

Implicit Function Theorem: Consider $\bm{\theta}_{J}^{*}(\bm{\theta}_{M})$
defined in Equation 11, with first-order optimality condition
$\frac{\partial\mathcal{L}_{so}(\bm{\theta}_{J}^{*},\bm{\theta}_{M})}{\partial\bm{\theta}_{J}}=0$,
$\begin{aligned}
&\frac{\partial}{\partial\bm{\theta}_{M}}\\!\left[\frac{\partial\mathcal{L}_{so}(\bm{\theta}_{J}^{*}(\bm{\theta}_{M}),\bm{\theta}_{M})}{\partial\bm{\theta}_{J}}\right]\\!=\\!0,\
\Rightarrow\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}\\!+\\!\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\frac{\partial(\bm{\theta}_{J}^{*}(\bm{\theta}_{M}))}{\partial\bm{\theta}_{M}}=0,\\\
&\Rightarrow\text{best-response Jacobian:
}\frac{\partial(\bm{\theta}_{J}^{*}(\bm{\theta}_{M}))}{\partial\bm{\theta}_{M}}=-\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{-1}\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}.\end{aligned}$
With Equation 12, we have hypergradient formulated by:

(14)
$\nabla_{\bm{\theta}_{M}}\mathcal{L}_{mo}=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{M}}-\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{-1}\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}.$

However, the inverse Hessian
$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{-1}$
in Equation 14 is hard to calculate. In Equation 13, BiSMO-FD employs finite
difference to naively approximate the inverse
$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{-1}=\xi\bm{\mathcal{I}}$,
where $\bm{\mathcal{I}}$ denotes the identity matrix. For a more precise
approximation of the inverse, we introduce two IFT-based methods: Neumann
series (BiSMO-NMN) and conjugate gradient (BiSMO-CG), to reformulate the
hypergradient.

####  3.2.2. BiSMO-NMN

######  Lemma 2.

Neumann series (lorraine2020optimizing, ): With a matrix $\bm{A}$ that
$\left\|\bm{\mathcal{I}}-\bm{A}\right\|<1$, we have
$\bm{A}^{-1}=\sum_{k=0}^{\infty}(\bm{\mathcal{I}}-\bm{A})^{k}$.

Based on 2, with small enough learning rate, the hypergradient in Equation 14
for BiSMO-NMN is formulated by:

(15) $\displaystyle\nabla_{\bm{\theta}_{M}}\mathcal{L}_{mo}$
$\displaystyle=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{M}}-\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\sum_{k=0}^{\infty}\left[\bm{\mathcal{I}}-\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{k}\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}.$

The approximation of the hypergradient
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}_{mo}$ can be derived by
considering only the first $K$ terms of the Neumann series, thereby avoiding
the need to calculate the inverse of the Hessian as:

(16)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{M}}-\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\sum_{k=0}^{K}\left[\bm{\mathcal{I}}-\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{k}\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}.$

####  3.2.3. BiSMO-CG:

Instead of calculating the Neumann series, another efficient way to
approximate the inverse Hessian is to solve the linear systems. Specifically,
$\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{-1}$
can be computed as the solution to the linear system
$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]\bm{w}=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}.$
The vector $\bm{w}$ can be obtained by solving the optimization problem:

(17)
$\min_{\bm{w}}\bm{w}^{\top}\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]\bm{w}-\bm{w}^{\top}\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}.$

The conjugate gradient (CG) algorithm is well-suited for this task, given its
efficient iteration complexity and use of Hessian-vector products (HVP) for
$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]\bm{w}$.
Such HVP can be obtained cheaply without explicitly forming or storing the
Hessian. The hypergradient in Equation 14 for BiSMO-CG is then computed as:

(18)

$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{CG}}_{mo}\\!=\\!\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{M}}\\!-\\!\left[\underset{\bm{w}}{\operatorname{argmin}}\left(\bm{w}^{\top}\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\bm{w}\\!-\\!\bm{w}^{\top}\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\right)\right]\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}.$

####  3.2.4. BiSMO-FD vs. BiSMO-NMN vs. BiSMO-CG vs. AM-SMO:

The optimization flow of BiSMO is demonstrated in Figure 2(b) and Algorithm 2.
When $K=0$, the
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}$ in Equation
16 reduces to match $\nabla_{\bm{\theta}_{M}}\mathcal{L}_{mo}^{\texttt{FD}}$
in Equation 13. In Algorithm 2, set $K=0;T=1;$ the FD can be executed through
the same process as NMN. Both NMN and CG use $T$ unroll SO steps to
approximate $\bm{\theta}_{J}^{*}$ (5), and employ HVP or Jacobian-vector
product (JVP) for computational acceleration. The key difference between them
lies in approximating the inverse Hessian: NMN uses the first $K$ terms of the
Neumann series (9), while CG applies $K$ CG steps (13). Compared to AM-SMO,
the unroll strategy in 5, due to the small T, substantially reduces runtime by
avoiding full SO cycle convergence. Furthermore, in 12 and 15, IFT-based
gradient fusion provides BiSMO with a global perspective and a more thorough
exploration of the solution space, thereby facilitating enhanced and
accelerated convergence.

Algorithm 2 Bilevel SMO via BiSMO-NMN, BiSMO-CG

1: Unroll step $T$, stepsizes $\xi_{J},\xi_{M}$, initializations
$\bm{\theta}_{J}^{0},\bm{\theta}_{M}^{0},\bm{w}_{0}$, term $K$

2: $\bm{\theta}_{J},\bm{\theta}_{M}$

3: $\bm{\theta}_{J},\bm{\theta}_{M}$ $\bm{\theta}_{J},\bm{\theta}_{M}$ while
not converged do $\bm{\theta}_{J},\bm{\theta}_{M}$
$\bm{\theta}_{J},\bm{\theta}_{M}$ while not converged do

4: $\bm{\theta}_{J},\bm{\theta}_{M}$ $\bm{\theta}_{J},\bm{\theta}_{M}$ while
not converged do $\bm{\theta}_{J},\bm{\theta}_{M}$
$\bm{\theta}_{J},\bm{\theta}_{M}$ while not converged do
$\bm{\theta}_{J},\bm{\theta}_{M}$ $\bm{\theta}_{J},\bm{\theta}_{M}$ while not
converged do for Inner step $t=1,\ldots,T$ do $\rhd$ Unroll $T$ steps of
inner-SO. $\bm{\theta}_{J},\bm{\theta}_{M}$ $\bm{\theta}_{J},\bm{\theta}_{M}$
while not converged do $\bm{\theta}_{J},\bm{\theta}_{M}$
$\bm{\theta}_{J},\bm{\theta}_{M}$ while not converged do for Inner step
$t=1,\ldots,T$ do $\rhd$ Unroll $T$ steps of inner-SO.

5: $\bm{\theta}_{J},\bm{\theta}_{M}$ $\bm{\theta}_{J},\bm{\theta}_{M}$ while
not converged do $\bm{\theta}_{J},\bm{\theta}_{M}$
$\bm{\theta}_{J},\bm{\theta}_{M}$ while not converged do for Inner step
$t=1,\ldots,T$ do $\rhd$ Unroll $T$ steps of inner-SO. Update
$\bm{\theta}_{J}^{t}\leftarrow\bm{\theta}_{J}^{t-1}-\xi_{J}\nabla_{\bm{\theta}_{J}}\mathcal{L}_{so}(\bm{\theta}_{J}^{t-1},\bm{\theta}_{M})$;
// Or Adam. Update
$\bm{\theta}_{J}^{t}\leftarrow\bm{\theta}_{J}^{t-1}-\xi_{J}\nabla_{\bm{\theta}_{J}}\mathcal{L}_{so}(\bm{\theta}_{J}^{t-1},\bm{\theta}_{M})$;
// Or Adam. Update
$\bm{\theta}_{J}^{t}\leftarrow\bm{\theta}_{J}^{t-1}-\xi_{J}\nabla_{\bm{\theta}_{J}}\mathcal{L}_{so}(\bm{\theta}_{J}^{t-1},\bm{\theta}_{M})$;
// Or Adam.

6: Update
$\bm{\theta}_{J}^{t}\leftarrow\bm{\theta}_{J}^{t-1}-\xi_{J}\nabla_{\bm{\theta}_{J}}\mathcal{L}_{so}(\bm{\theta}_{J}^{t-1},\bm{\theta}_{M})$;
// Or Adam. Approximate $\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$;
Re-initialize $\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$; Approximate
$\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$; Re-initialize
$\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$;

7: Approximate $\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$; Re-
initialize $\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$; Approximate
$\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$; Re-initialize
$\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$; Approximate
$\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$; Re-initialize
$\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$; if BiSMO-NMN then $\rhd$
Hypergradient via BiSMO-NMN. Approximate
$\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$; Re-initialize
$\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$; Approximate
$\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$; Re-initialize
$\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$; if BiSMO-NMN then $\rhd$
Hypergradient via BiSMO-NMN.

8: Approximate $\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$; Re-
initialize $\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$; Approximate
$\bm{\theta}_{J}^{*}\leftarrow\bm{\theta}_{J}^{t}$; Re-initialize
$\bm{\theta}_{J}^{0}\leftarrow\bm{\theta}_{J}^{t}$; if BiSMO-NMN then $\rhd$
Hypergradient via BiSMO-NMN. 1) Get $K$ Neumann series via HVP:
$\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\sum_{k=0}^{K}\left[\bm{\mathcal{I}}-\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{k}$;
1) Get $K$ Neumann series via HVP:
$\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\sum_{k=0}^{K}\left[\bm{\mathcal{I}}-\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{k}$;

9: 1) Get $K$ Neumann series via HVP:
$\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}}\sum_{k=0}^{K}\left[\bm{\mathcal{I}}-\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]^{k}$;
2) Get Jacobian-vector product JVP in Equation 16; 2) Get Jacobian-vector
product JVP in Equation 16;

10: 2) Get Jacobian-vector product JVP in Equation 16; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16;

11: 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; if BiSMO-CG then $\rhd$ Hypergradient via BiSMO-CG. 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; if BiSMO-CG then $\rhd$ Hypergradient via BiSMO-CG.

12: 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{NMN}}_{mo}\leftarrow$
Equation 16; if BiSMO-CG then $\rhd$ Hypergradient via BiSMO-CG. 1) Solve
$\bm{w}_{K}$ from
$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]\bm{w}=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}},$
via $K$ steps of CG starting from $\bm{w}_{0}$; then re-initialize
$\bm{w}_{0}\leftarrow\bm{w}_{K}$; 1) Solve $\bm{w}_{K}$ from
$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]\bm{w}=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}},$
via $K$ steps of CG starting from $\bm{w}_{0}$; then re-initialize
$\bm{w}_{0}\leftarrow\bm{w}_{K}$;

13: 1) Solve $\bm{w}_{K}$ from
$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{J}\partial\bm{\theta}_{J}}\right]\bm{w}=\frac{\partial\mathcal{L}_{mo}}{\partial\bm{\theta}_{J}},$
via $K$ steps of CG starting from $\bm{w}_{0}$; then re-initialize
$\bm{w}_{0}\leftarrow\bm{w}_{K}$; 2) Get Jacobian-vector product
JVP:$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}\right]\bm{w}$;
2) Get Jacobian-vector product
JVP:$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}\right]\bm{w}$;

14: 2) Get Jacobian-vector product
JVP:$\left[\frac{\partial^{2}\mathcal{L}_{so}}{\partial\bm{\theta}_{M}\partial\bm{\theta}_{J}}\right]\bm{w}$;
3) $\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{CG}}_{mo}\leftarrow$
Equation 18; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{CG}}_{mo}\leftarrow$
Equation 18; 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{CG}}_{mo}\leftarrow$
Equation 18;

15: 3)
$\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}^{\texttt{CG}}_{mo}\leftarrow$
Equation 18; Update
$\bm{\theta}_{M}\leftarrow\bm{\theta}_{M}-\xi_{M}\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}_{mo}$;
// Or Adam. Update
$\bm{\theta}_{M}\leftarrow\bm{\theta}_{M}-\xi_{M}\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}_{mo}$;
// Or Adam. Update
$\bm{\theta}_{M}\leftarrow\bm{\theta}_{M}-\xi_{M}\nabla_{\bm{\theta}_{M}}\tilde{\mathcal{L}}_{mo}$;
// Or Adam.

##  4\. Experiments

The BiSMO is implemented in PyTorch framework and tested on an Nvidia RTX4090
GPU card across three datasets as listed in Table 2. The hyperparameters
settings are as follows
$\gamma\\!=\\!1000;\eta\\!=\\!3000;\lambda\\!=\\!193nm;N\\!A\\!=\\!1.35;\sigma_{o}\\!=\\!0.95;\sigma_{i}\\!=\\!0.63.Q\\!=\\!24;N_{j}\\!=\\!35;N_{m}\\!=\\!2048,\alpha_{m}\\!=\\!9;m_{0}\\!=\\!1;\alpha_{j}\\!=\\!2;j_{0}\\!=\\!5;\beta\\!=\\!30;P\\!=\\!256;\xi\\!=\\!\xi_{M}\\!=\\!\xi_{J}\\!=\\!0.1;K\\!=\\!5;T\\!=\\!3$.
All the tiles are converted to $2048\times 2048$-pixel images.

Table 2. Details of the Dataset.

Dataset | From |  Area† | Test num. | Layer |  CD‡ | tile
---|---|---|---|---|---|---
ICCAD13 | (OPC-DAC2014-Gao, ; OPC-DAC2018-Yang, ) | 202655 | 10 | Metal |  32$nm$ |  4$\mu m^{2}$
ICCAD-L | (OPC-ICCAD2020-NeuralILT, ; multiilt2023sun, ) | 475571 | 10 | Metal |  32$nm$ |  4$\mu m^{2}$
ISPD19 | (ispd2019-benchmark, ) | 698743 | 100 | Metal+Via |  28$nm$ |  4$\mu m^{2}$
Area†: average area: unit $nm^{2}$; CD‡: critical dimension. |

###  4.1. Results Comparison with SOTA

Table 3. Result comparison with SOTA.

Bench | MO | AM-SMO |  BiSMO (Ours)
---|---|---|---
NILT (OPC-ICCAD2020-NeuralILT, ) | DAC23-MILT (multiilt2023sun, ) | Abbe-MO(Ours) | Abbe-Hopkins⋆ (ding2020gradient, ) | Abbe-Abbe (sampleICC2022sun, ) | BiSMO-FD | BiSMO-CG | BiSMO-NMN
L2 | PVB | L2 | PVB | L2 | PVB | L2 | PVB | L2 | PVB | L2 | PVB | L2 | PVB | L2 | PVB
ICCAD13 | 37515 | 50964 | 28362 | 40044 | 20419 | 29697 | 27299 | 37278 | 17539 | 23944 | 13828 | 17872 | 13603 | 16274 | 13059 | 15839
ICCAD-L | 71570 | 108162 | 53143 | 87010 | 44478 | 66092 | 48879 | 77062 | 40455 | 58560 | 29779 | 42643 | 29762 | 40543 | 28946 | 38706
ISPD19 | 97891 | 119732 | 85234 | 105592 | 61374 | 93132 | 79634 | 97073 | 55588 | 84402 | 39959 | 64211 | 39488 | 61190 | 38737 | 59832
Average | 68992 | 92953 | 55580 | 77549 | 42090 | 62974 | 51937 | 70471 | 37861 | 55635 | 27855 | 41576 | 27618 | 39336 | 26914 | 38126
Ratio | 2.56 | 2.44 | 2.07 | 2.03 | 1.56 | 1.65 | 1.93 | 1.85 | 1.41 | 1.46 | 1.03 | 1.09 | 1.03 | 1.03 | 1.00 | 1.00
Abbe-Hopkins⋆ (ding2020gradient, ): AM-SMO employs Abbe model for SO and Hopkins model for MO. L2 and PVB unit: $nm^{2}$. |

Table 4. EPE and runtime comparison.

| MO | AM-SMO |  | BiSMO |
---|---|---|---|---|---
NILT (OPC-ICCAD2020-NeuralILT, ) |  DAC23 (multiilt2023sun, ) |  Abbe -MO |  A$\sim$H⋆ (ding2020gradient, ) |  A$\sim$A$\dagger$ (sampleICC2022sun, ) | FD | CG | NMN
EPE avg. | 10.1 | 3.6 | 2.8 | 3.3 | 2.4 | 1.8 | 1.6 | 1.6
Ratio | 6.2 | 2.2 | 1.7 | 2.0 | 1.5 | 1.1 | 1.0 | 1.0
TAT‡ avg. | 12.4 | 3.8 | 11.7 | 287 | 122.5 | 12.6 | 15.3 | 14.7
Ratio | 0.84 | 0.26 | 0.80 | 19.52 | 8.33 | 0.86 | 1.04 | 1.00
A$\sim$H⋆: Abbe-Hopkins; A$\sim$A†: Abbe-Abbe; TAT‡: Turn around time (s).

We have conducted a comparative analysis of the performance between our Abbe-
based MO and the previous SOTA MO methods DAC23-MILT (multiilt2023sun, ) and
NILT (OPC-ICCAD2020-NeuralILT, ) in Table 3 and Table 4. Furthermore, we
compare the performance of BiSMO with the previous SOTA AM-SMO
(sampleICC2022sun, ; ding2020gradient, ). AM-SMO is implemented in two ways:
one involves hybrid Abbe-SO and Hopkins-MO (ding2020gradient, ), while the
other employs the Abbe model for both SO and MO (sampleICC2022sun, ).

To highlight advantages of BiSMO, Figure 3 shows log-scaled loss
$\log(\mathcal{L}_{smo})$ convergence compared to SOTA MO (multiilt2023sun, )
and AM-SMO (sampleICC2022sun, ) using random cases from test datasets in Table
2, with a 0.01 learning rate. Result samples are depicted in Figure 4,
appropriately scaled and cropped to enhance visualization.

---

$0$$20$$80$$100$$0.2$$0.7$step$\log(\mathcal{L}_{smo})$ ((a))

$0$$20$$80$$100$0.10.6step$\log(\mathcal{L}_{smo})$ ((b))

$0$$20$$80$$100$0.10.5step$\log(\mathcal{L}_{smo})$ ((c))

$0$$20$$80$$100$0.10.5step$\log(\mathcal{L}_{smo})$DAC23(multiilt2023sun,
)Abbe-MOAM-SMO(sampleICC2022sun, )BiSMO-FDBiSMO-CGBiSMO-NMN ((d))

Figure 3. Loss comparison between different MO methods (dashed lines) and SMO
methods (solid lines).

Effectiveness of Abbe-MO: Table 3 and Table 4 demonstrate that our Abbe-MO
achieved a 25% reduction in L2, 19% in PVB, and 24% in EPE when compared to
the SOTA MO DAC23-MILT (multiilt2023sun, ). Furthermore, within the AM-SMO,
the Abbe-based SMO(sampleICC2022sun, ) surpassed the Abbe-Hopkins hybrid
SMO(ding2020gradient, ) by decreasing L2 by 28%, PVB by 21%, and EPE by 29%.
In Figure 3, all cases indicate that Abbe-MO converges more rapidly and
effectively than SOTA MO (multiilt2023sun, ). This can be attributed to the
fact that truncated decomposition in Hopkins’ approach leads to a loss of
accuracy in lithography, thereby allowing our lossless Abbe method to achieve
superior MO and SMO results.

BiSMO vs. AM-SMO vs. MO: Table 3, Table 4 and Figure 3 reveal that BiSMO
variants significantly outperform AM-SMO(ding2020gradient, ; sampleICC2022sun,
) in error reduction, with BiSMO-NMN achieving decreases of 41% in L2, 46% in
PVB, and 37% in EPE, and even the basic BiSMO-FD showing reductions of 36% in
L2, 34% in PVB, and 27% in EPE. Figure 3 illustrates AM-SMO’s
(sampleICC2022sun, ) ‘zigzag’ loss curve, a result of its alternating
optimization, which ultimately settles below Abbe-MO yet above BiSMO variants.
This suggests that while AM-SMO’s broader solution space improves outcomes
compared to pure MO. However, its alternating approach risks entrapment in
local minima, hindering it from reaching the optimal results achievable by
BiSMO. This fact underscores the superior performance of the BiSMO method.
Additionally, BiSMO demonstrates a significant improvement over the SOTA MO
(multiilt2023sun, ), achieving a $\sim$50% reduction in all error metrics.

Runtime comparison: In our implementation, Abbe-MO and Hopkins-MO have been
accelerated to average 0.16s and 0.12s per MO iteration, aligning with
theoretical derivations in Table 1. As shown in Table 4, for parity, we have
applied GPU acceleration to AM-SMO (ding2020gradient, ; sampleICC2022sun, )
with settings identical to BiSMO. In their original implementations, the
runtime for (sampleICC2022sun, ) was 910s, and (ding2020gradient, ) was 69
minutes. Utilizing our accelerated Abbe imaging, this has been accelerated to
122.5s and 287s, respectively. Despite these improvements, BiSMO, leveraging
hypergradient, still achieves faster convergence than (ding2020gradient, ;
sampleICC2022sun, ), boosting throughput by 8.3 times compared to the Abbe-
based AM-SMO (sampleICC2022sun, ). Furthermore, it’s about 19.5 times quicker
than the Abbe-Hopkins hybrid AM-SMO (ding2020gradient, ), which is slowed down
by its complex iterative TCC generation and decomposition.

Figure 4. Result samples from ICCAD13 and ISPD19 datasets.

###  4.2. Ablation Study

BiSMO-FD vs. BiSMO-NMN vs. BiSMO-CG: Figure 3, Table 3, and Table 4 clearly
show that NMN typically outperforms other methods, followed by CG, with FD
being relatively weaker among all BiSMO variants. The relative instability of
CG is indicated by its largest standard deviation (STD) in Figure 5.
Meanwhile, FD boasts the shortest runtime (Table 4), and CG’s advantage lies
in outperforming NMN in some cases, as shown in Figure 3(d).

---

20305060125step$\mathcal{L}_{smo}$ ((a))

20305060367step$\mathcal{L}_{smo}$BiSMO-FD meanBiSMO-CG meanBiSMO-NMN mean
((b))

Figure 5. Mean and STD of (a) ICCAD (b) ICCAD-L datasets.

##  5\. Conclusion

In this paper, we first establish a unified SMO framework utilizing Abbe-
imaging, enabling simultaneous SO and MO gradient computation with enhanced
lithographic precision and rapid calculation. Building on this foundation,
BiSMO is introduced, conceptualizing SMO as a bilevel problem and proposing
three innovative methods for calculating source-mask best-response gradients,
effectively addressing bilevel SMO challenges. BiSMO’s gradient-based
approach, with its global perspective and improved exploration of solution
space facilitates navigation out of local minima, ensuring better and faster
converging SMO outcomes. This method surpasses traditional AM-SMO limitations,
positioning bilevel SMO as a future promising mainstream approach in the
field.

## References

  * (1) C. H. Wallace, P. A. Nyhus, and S. S. Sivakumar, “Sub-resolution assist features,” Dec. 15 2009, US Patent.
  * (2) J.-R. Gao, X. Xu, B. Yu, and D. Z. Pan, “MOSAIC: Mask optimizing solution with process window aware inverse correction,” in _Proc. DAC_ , 2014.
  * (3) P. Evanschitzky, A. Erdmann, and T. Fuehner, “Extended abbe approach for fast and accurate lithography imaging simulations,” in _25th European Mask and Lithography Conference_ , 2009.
  * (4) N. Cobb, “Sum of coherent system decomposition by SVD,” _Berkeley CA_ , 1995\.
  * (5) H. Yang, S. Li, Y. Ma, B. Yu, and E. F. Young, “GAN-OPC: Mask optimization with lithography-guided generative adversarial nets,” in _Proc. DAC_ , 2018\.
  * (6) G. Chen, W. Chen, Y. Ma, H. Yang, and B. Yu, “DAMO: Deep agile mask optimization for full chip scale,” in _Proc. ICCAD_ , 2020.
  * (7) B. Jiang, L. Liu, Y. Ma, H. Zhang, E. F. Y. Young, and B. Yu, “Neural-ILT: Migrating ILT to nerual networks for mask printability and complexity co-optimizaton”,” in _Proc. ICCAD_ , 2020.
  * (8) G. Chen, Z. Yu, H. Liu, Y. Ma, and B. Yu, “DevelSet: Deep neural level set for instant mask optimization,” in _Proc. ICCAD_ , 2021.
  * (9) Z. Yu, G. Chen, Y. Ma, and B. Yu, “A gpu-enabled level set method for mask optimization,” in _Proc. DATE_ , 2021.
  * (10) S. Sun, F. Yang, B. Yu, L. Shang, and X. Zeng, “Efficient ILT via multi-level lithography simulation,” in _Proc. DAC_ , 2023.
  * (11) Z. Wang, X. Ma, R. Chen, S. Zhang, and G. R. Arce, “Fast pixelated lithographic source and mask joint optimization based on compressive sensing,” _IEEE Transactions on Computational Imaging_ , vol. 6, pp. 981–992, 2020.
  * (12) Y. Sun, Y. Li, G. Liao, M. Yuan, P. Wei, Y. Li, L. Zou, and L. Liu, “Sampling-based imaging model for fast source and mask optimization in immersion lithography,” _Appl. Opt._ , 2022.
  * (13) M. Ding, Z. Niu, F. Zhang, L. Zhu, W. Shi, A. Zeng, and H. Huang, “Gradient-based source mask and polarization optimization with the hybrid hopkins–abbe model,” _Journal of Micro/Nanolithography, MEMS, and MOEMS_ , 2020.
  * (14) J. Lorraine, P. Vicol, and D. Duvenaud, “Optimizing millions of hyperparameters by implicit differentiation,” in _International Conference on Artificial Intelligence and Statistics_ , 2020.
  * (15) M. Zhang, S. W. Su, S. Pan, X. Chang, E. M. Abbasnejad, and R. Haffari, “iDARTS: Differentiable architecture search with stochastic implicit gradients,” in _Proc. ICML_ , 2021.
  * (16) S. K. Choe, S. V. Mehta, H. Ahn, W. Neiswanger, P. Xie, E. Strubell, and E. Xing, “Making scalable meta learning practical,” in _Proc. NeurIPS_ , 2023.
  * (17) S. Banerjee, Z. Li, and S. R. Nassif, “ICCAD-2013 CAD contest in mask optimization and benchmark suite,” in _Proc. ICCAD_ , 2013.
  * (18) W.-H. Liu, S. Mantik, W.-K. Chow, Y. Ding, A. Farshidi, and G. Posser, “ISPD 2019 initial detailed routing contest and benchmark with advanced routing rules,” in _Proc. ISPD_ , 2019.


# Impact of high-rank excitations on accuracy of the unitary coupled cluster
downfolding formalism

Karol Kowalski<EMAIL_ADDRESS>Physical Sciences Division, Pacific
Northwest National Laboratory, Richland, Washington 99354, United States  Bo
Peng  Physical Sciences Division, Pacific Northwest National Laboratory,
Richland, Washington 99354, United States  Nicholas P. Bauman  Physical
Sciences Division, Pacific Northwest National Laboratory, Richland, Washington
99354, United States

###### Abstract

In this paper, we evaluate the accuracy of the Hermitian form of the
downfolding procedure utilizing the double unitary coupled cluster Ansatz
(DUCC) on the H6 and H8 benchmark systems. The computational infrastructure
employs the occupation-number-representation codes to construct the matrix
representation of arbitrary second-quantized operators, enabling the exact
representation of exponentials of various operators. The tests utilize
external excitations estimated from standard single-reference coupled cluster
methods (SR-CC) to demonstrate that higher-rank SR-CC external amplitudes were
necessary to describe the energies in the strongly correlated regime
adequately. We show that this approach can offset problems of the
corresponding SR-CC theories associated with losing the variational character
of corresponding energies.

##  I Introduction

Applying many-body methods for dimensionality/cost reduction (DCR) of ab-
initio formulations is imperative in expanding the range of system sizes
amenable to accurate many-body formulations in chemistry and material
sciences. Additionally, these techniques are vital in effectively using early
quantum computing resources, commonly referred to as the noisy intermediate-
scale quantum devices (NISQ). Takeshita  _et al._ (2020); Şahinoğlu and Somma
(2021); Liu, Li, and Yang (2022); Huang, Li, and Evangelista (2023) DCR
methods primarily focus on minimizing the number of qubits required to
represent a given quantum problem. One should mention several techniques
developed to take full advantage of the ubiquitous Variational Quantum
Eigensolvers (VQE) approach Peruzzo  _et al._ (2014); McClean  _et al._
(2016); Romero  _et al._ (2018); Kandala  _et al._ (2017, 2019); Izmaylov  _et
al._ (2019); Lang, Ryabinkin, and Izmaylov (2021); Grimsley  _et al._ (2019a,
b); McArdle  _et al._ (2020); Kirby and Love (2021); Tilly  _et al._ (2022) in
addressing problems beyond the situation where few electrons are correlated.

In the context of quantum algorithms for quantum chemistry, the utilization of
DCR techniques is linked to the partitioning of electron correlation effects
into static and dynamic partitions. In terms of methodology, the coupled
cluster (CC) formalism Coester (1958); Coester and Kummel (1960); Čížek
(1966); Paldus, Čížek, and Shavitt (1972); Mukherjee, Moitra, and Mukhopadhyay
(1975); Adams, Jankowski, and Paldus (1979); Purvis and Bartlett (1982);
Jeziorski and Monkhorst (1981); Lyakh  _et al._ (2012) offers an effective
means of describing these effects in a many-body language. Although static
effects can be incorporated for small-scale systems using presently accessible
quantum hardware, the implementation of necessary dynamical correlation
effects, which typically involve numerous parameters with minute values,
remains beyond the scope of contemporary quantum computing technologies.

We recently introduced and tested downfolding techniques based on the double
unitary coupled cluster Ansatz (DUCC) Bauman  _et al._ (2019) to address the
abovementioned problem. The downfolding procedure utilizes the properties of
the ground-state DUCC Ansatz, which in analogy to single-reference sub-system
embedding sub-algebras (SES-CC), Kowalski (2018, 2021) allows to construct of
effective Hamiltonians that integrate out-of-active-space degrees of freedom
usually identified with dynamical amplitudes. In contrast to the SES-CC
approach, the DUCC formalism yields the Hermitian form of the effective
Hamiltonian in the active space.

The DUCC-driven downfolded Hamiltonians are critical components of hybrid
computing. Classical computing resources are employed to calculate the second
quantized form of effective Hamiltonians, and quantum computing is invoked to
diagonalize them in active spaces that ideally match the available quantum
resources. This type of approach provides much needed algorithmic transition
mechanism from current NISQ technologies Preskill (2018) to mature error-
corrected quantum computers of the future, where the size of the active space
is adjusted to the available quantum resources. For this purpose, several
approximations were tested to validate the efficiency of the downfolding
procedure. These approximations, due to the non-commutativity of the
components defining DUCC cluster operators, were based on the finite low-rank
commutator expansions, the limited rank of interactions included in the
downfolded Hamiltonians (one- and two-body interactions), and simple form of
the external amplitudes extracted from the single-reference CC (SR-CC) model
with singles and doubles (CCSD).Purvis and Bartlett (1982)

Our team has recently developed a novel full configuration interaction (FCI)
code called stringMB, which employs a string-based approach to emulate quantum
systems and represent operators in matrix form. This code has been integrated
into the NWChem software, enabling us to (1) work with the exact
representations of operator exponents and (2) leverage various sources for
external CC amplitudes. Consequently, we have a unique opportunity to study
the exact nature of downfolded Hamiltonians using the DUCC method. In this
study, we investigate the impact of higher-rank external excitations obtained
through CCSD,Purvis and Bartlett (1982) CCSDT,Scuseria and Schaefer (1988);
Noga and Bartlett (1987, 1988) and CCSDTQ Oliphant and Adamowicz (1991a);
Piecuch and Adamowicz (1994); Kucharski and Bartlett (1991) simulations, as
well as the active space size, on the accuracy of ground-state energies for
small benchmark systems H6 and H8 representing linear chains of hydrogen
atoms.

##  II Theory

The DUCC formulations have been amply discussed in recent papers (see
Refs.Bauman  _et al._ (2019); Kowalski and Bauman (2020); Bauman and Kowalski
(2022a)). Here we overview only the salient features of these approaches.
While the SES-CC technique Kowalski (2018, 2021) forms the basis for non-
Hermitian downfolding, the DUCC expansions provide its Hermitian formulations.
The Hermitian form of the downfolded Hamiltonian is obtained as a consequence
of utilizing active-space-dependent DUCC representation of the wave function

$|\Psi\rangle=e^{\sigma_{\rm ext}}e^{\sigma_{\rm int}}|\Phi\rangle\;,$ (1)

where $\sigma_{\rm ext}$ and $\sigma_{\rm int}$, referred to as the external
and internal cluster operators, are general-type anti-Hermitian operators

$\displaystyle\sigma_{\rm int}^{\dagger}$ $\displaystyle=$
$\displaystyle-\sigma_{\rm int}\;,$ (2) $\displaystyle\sigma_{\rm
ext}^{\dagger}$ $\displaystyle=$ $\displaystyle-\sigma_{\rm ext}\;.$ (3)

In analogy to the non-Hermitian case, the $\sigma_{\rm ext}$ and $\sigma_{\rm
int}$ operators are defined by parameters carrying active spin-orbital labels
only and at least one in-active spin-orbital label, respectively. The DUCC
Ansatz falls into a broad class of active space coupled cluster methods.
Oliphant and Adamowicz (1991b, 1992); Piecuch, Oliphant, and Adamowicz (1993)

The use of the DUCC Ansatz (1), in analogy to the SES-CC case, leads to an
alternative way of determining energy, which can be obtained by solving
active-space Hermitian eigenvalue problem:

$H^{\rm eff}e^{\sigma_{\rm int}}|\Phi\rangle=Ee^{\sigma_{\rm
int}}|\Phi\rangle,$ (4)

where

$H^{\rm eff}=(P+Q_{\rm int})\bar{H}_{\rm ext}(P+Q_{\rm int})$ (5)

and

$\bar{H}_{\rm ext}=e^{-\sigma_{\rm ext}}He^{\sigma_{\rm ext}}.$ (6)

When the external cluster amplitudes are known (or can be effectively
approximated), the energy (or its approximation) can be calculated by
diagonalizing the Hermitian effective/downfolded Hamiltonian (5) in the active
space using various quantum or classical diagonalizers. The $Q_{\rm int}$
operator is a projection onto excited (with respect to $|\Phi\rangle$) excited
configurations in complete active space (CAS) and the projection onto the
reference function is denoted as $P$.

For quantum computing applications second-quantized representation of $H^{\rm
eff}$ is required. In the light of the non-commuting character of components
defining $\sigma_{\rm ext}$ operator, to this end, one has to rely on the
finite-rank commutator expansions, i.e.,

$\bar{H}_{\rm ext}\simeq H+\sum_{i=1}^{\rm
Max_{R}}\frac{1}{i!}[\ldots[H,\sigma_{\rm ext}],\ldots],\sigma_{\rm
ext}]_{i}\;,$ (7)

where ${\rm Max_{R}}$ stands for the length of commutator expansion. Due to
the numerical costs associated with the contractions of multi-dimensional
tensors, only approximations based on including low-rank commutators are
feasible. In recent studies, approximations based on single, double, and part
of triple commutators were explored where one- and two-body interactions were
retained in the second quantized form of $H^{\rm eff}$. In practical
applications, one also has to determine the approximate form of $\sigma_{\rm
ext}$. For practical reasons we used the following approximation

$\sigma_{\rm ext}\simeq T_{\rm ext}-T_{\rm ext}^{\dagger}\;,$ (8)

where $T_{\rm ext}$ can be defined through the external parts of the typical
SR-CC cluster operators.

Given the progress achieved in the development of the Hermitian form of the
downfolded Hamiltonians, addressing two pressing questions (1) what is the
impact of the choice of the $T_{\rm ext}$ on the quality of ground-state
energy of $H^{\rm eff}$? and (2) what are the energy values corresponding to
the untruncated (exact) form of the $H^{\rm eff}$? play a pivotal role in
further understanding and advances of CC downfolding techniques. We answer
these questions using stringMB code that allows us to deal with the exact
matrix representations of second quantized operators and their functions in
the FCI space.

##  III Implementation

For interacting fermionic systems,the action of the creation/annihilation
operators for the electron in $p$-th spin-orbital ($a_{p}/a_{p}^{\dagger}$) on
the Slater determinants can be conveniently described using occupation number
representation, where each Slater determinant is represented as a vector

$|n_{M}\;n_{M-1}\;\ldots\;n_{i+1}\;n_{i}\;n_{i-1}\;\ldots\;n_{1}\rangle$ (9)

where occupation numbers $n_{i}$ are equal to either 1 (electron occupies
$i$-th spin orbital) or 0 (no electron is occupying $i$-th spin orbital). In
(9), $M$ stands for the total number of spin-orbitals used to describe quantum
system and $M=2n$, where $n$ is the number of orbitals.

The following formulas give the non-trivial action of creation/annihilation
operators on the state vectors

$\displaystyle
a_{i}^{\dagger}|n_{M}\;n_{M-1}\;\ldots\;n_{i+1}\;0\;n_{i-1}\;\ldots
n_{1}\;\rangle$ $\displaystyle=$
$\displaystyle(-1)^{\sum_{k=1}^{i-1}n_{k}}|n_{M}\;n_{M-1}\;\ldots\;n_{i+1}\;1\;n_{i-1}\;\ldots
n_{1}\;\rangle\;\;\;\;\;\;$ (10) $\displaystyle
a_{i}|n_{M}\;n_{M-1}\;\ldots\;n_{i+1}\;1\;n_{i-1}\;\ldots n_{1}\;\rangle$
$\displaystyle=$
$\displaystyle(-1)^{\sum_{k=1}^{i-1}n_{k}}|n_{M}\;n_{M-1}\;\ldots\;n_{i+1}\;0\;n_{i-1}\;\ldots
n_{1}\;\rangle.$ (11)

Using the occupation-number representation, the stringMB code allows one to
construct a matrix representation (${\bf A}$) of general second-quantized
operators $A$, where $A$ can be identified with electronic Hamiltonian, the
external part of the cluster operator $T_{\rm ext}$, and exponents of $T_{\rm
ext}-T_{\rm ext}^{\dagger}$, i.e.,

$\displaystyle H$ $\displaystyle\rightarrow$ $\displaystyle{\bf H}\;\;,$ (12)
$\displaystyle T_{\rm ext}$ $\displaystyle\rightarrow$ $\displaystyle{\bf
T}_{\rm ext}\;\;,$ (13) $\displaystyle e^{\sigma_{\rm ext}}\simeq e^{T_{\rm
ext}-T_{\rm ext}^{\dagger}}$ $\displaystyle\rightarrow$ $\displaystyle e^{{\bf
T}_{\rm ext}-{\bf T}_{\rm ext}^{\dagger}},$ (14) $\displaystyle
e^{-\sigma_{\rm ext}}\simeq e^{-(T_{\rm ext}-T_{\rm ext}^{\dagger})}$
$\displaystyle\rightarrow$ $\displaystyle e^{-({\bf T}_{\rm ext}-{\bf T}_{\rm
ext}^{\dagger})},$ (15) $\displaystyle\bar{H}_{\rm ext}$
$\displaystyle\rightarrow$ $\displaystyle\bar{\bf H}_{\rm ext},$ (16)
$\displaystyle H^{\rm eff}$ $\displaystyle\rightarrow$ $\displaystyle{\bf
H}^{\rm eff}\;.$ (17)

Moreover, the stringMB can extract sub-blocks of matrices or their products
corresponding to arbitrary active space. This feature is used to form matrix
representations of the effective Hamiltonians $H^{\rm eff}$.

##  IV Results

Owing to the memory requirements (associated with the storage of matrix
representations of the operator) of the stringMB code, we can deal with
relatively small systems yet epitomizing situations encountered in the
calculations for larger systems and processes. For this reason, we employed
ubiquitous models corresponding to the linear chains of hydrogen atoms: H6 and
H8 models. By varying neighboring hydrogens distance ($R_{\rm H-H}$), one can
smoothly transition from the single-reference character of the ground-state
wave function for smaller $R_{\rm H-H}$ distances ($\simeq$ 2.0 a.u. or less)
to quasi-degenerate regime (2.75 and 3.00 a.u.) - a typical situation
encountered in bond breaking/forming processes. These benchmark systems are
commonly used in accuracy studies of ab-into methodologies to deal with the
strong correlation effects. We used STO-3G basis set Hehre, Stewart, and Pople
(1969) and in all calculations, we used restricted Hartree-Fock (RHF)
molecular bases composed of 6 and 8 molecular orbitals for H6 and H8 systems,
respectively. For the strong correlation regime ($R_{\rm H-H}=3.0$ a.u.),
there is no obvious choice of the active space for the ground-state problem;
according to the wave function analysis for large separations, all orbitals
fall into the category of ”important” or ”active” orbitals. The results of our
simulations for H6 and H8 systems are summarized in Tables 1 \- 6.

In Table 1, we collected CC and DUCC results obtained for active space defined
by two highest occupied orbitals (orbitals 2 and 3) and two lowest virtual
orbitals (orbitals 4 and 5). The canonical CCSD, CCSDT, and CCSDTQ energies
(corresponding columns are denoted as SD, SDT, and SDTQ) are compared to the
FCI ones and lowest eigenvalues of the downfolded Hamiltonian in the active
space for CCSD (DUCC-SD), CCSDT (DUCC-SDT), and CCSDTQ (DUCC-SDTQ) sources of
the external amplitudes $T_{\rm ext}$ to calculate $\sigma_{\rm
ext}(\mathfrak{h})$ according to formula (8). While for weakly correlated
regimes, all choices of $T_{\rm ext}$ result in comparable accuracy with
respect to the FCI energies, for larger $R_{\rm H-H}$ distances the quality of
$T_{\rm ext}$ matters. For example, for larger distances, only $T_{\rm ext}$
defined at the CCSDTQ level yields DUCC-SDTQ energies of the CCSDTQ/FCI
quality for all geometries. Other approaches, CCSD and CCSDT, for larger
values of $R_{\rm H-H}$, provide energies significantly below the FCI ones.
However, it is interesting to notice that the corresponding eigenvalues of
DUCC-SD and DUCC-SDT can reinstate the variational character of ground-state
energy despite using $T_{\rm ext}$ stemming from CCSD and CCSDT calculations,
which is an important property of the DUCC approach.

$R_{\rm H-H}$ | FCI | SD | SDT | SDTQ | DUCC-SD | DUCC-SDT | DUCC-SDTQ
---|---|---|---|---|---|---|---
1.50 | -3.199566 | -3.199332 | -3.199601 | -3.199566 | -3.199324 | -3.199562 | -3.199566
1.75 | -3.245936 | -3.245603 | -3.246054 | -3.245936 | -3.245547 | -3.245923 | -3.245936
2.00 | -3.217699 | -3.217277 | -3.218047 | -3.217699 | -3.217040 | -3.217655 | -3.217697
2.25 | -3.156624 | -3.156266 | -3.157559 | -3.156621 | -3.155447 | -3.156484 | -3.156618
2.50 | -3.085398 | -3.085691 | -3.087713 | -3.085380 | -3.083217 | -3.084962 | -3.085374
2.75 | -3.016841 | -3.019512 | -3.022159 | -3.016770 | -3.012642 | -3.015537 | -3.016758
3.00 | -2.957646 | -2.967326 | -2.969163 | -2.957405 | -2.948732 | -2.953850 | -2.957384
Table 1: Comparison of energies of the downfolded Hamiltonians for the linear H6 system in the STO-3G basis set based on various sources of the external amplitudes $T_{\rm ext}$ used to approximate the $\sigma_{\rm ext}$ operator ($\sigma_{\rm ext}\simeq T_{\rm ext}-T_{\rm ext}^{\dagger}$). All simulations used restricted Hartree-Fock molecular orbitals 2,3, and 4,5 as active occupied and virtual orbitals, respectively. In the linear chain of the H atoms, the geometry is defined by the distance between neighboring hydrogen atoms ($R_{\rm H-H}$) in a.u. $R_{\rm H-H}$ | FCI | SD | SDT | SDTQ | DUCC-SD | DUCC-SDT | DUCC-SDTQ
---|---|---|---|---|---|---|---
1.50 | -4.235775 | -4.235111 | -4.235846 | -4.235775 | -4.235071 | -4.235757 | -4.235774
1.75 | -4.315273 | -4.314347 | -4.315504 | -4.315273 | -4.314173 | -4.315222 | -4.315271
2.00 | -4.286011 | -4.284844 | -4.286688 | -4.286013 | -4.284235 | -4.285862 | -4.286005
2.25 | -4.208339 | -4.207232 | -4.210169 | -4.208337 | -4.205334 | -4.207876 | -4.208316
2.50 | -4.114829 | -4.115000 | -4.119502 | -4.114795 | -4.109473 | -4.113350 | -4.114739
2.75 | -4.023783 | -4.029321 | -4.035510 | -4.023578 | -4.013082 | -4.018712 | -4.023447
3.00 | -3.944748 | -3.972672 | -3.978401 | -3.943920 | -3.912005 | -3.921323 | -3.943614
Table 2: Comparison of energies of the downfolded Hamiltonians for the linear
H8 system in the STO-3G basis set based on various sources of the external
amplitudes $T_{\rm ext}$ used to approximate the $\sigma_{\rm ext}$ operator
($\sigma_{\rm ext}\simeq T_{\rm ext}-T_{\rm ext}^{\dagger}$). All simulations
used restricted Hartree-Fock molecular orbitals 2,3, and 4,5 as active
occupied and virtual orbitals, respectively. In the linear chain of the H
atoms, the geometry is defined by the distance between neighboring hydrogen
atoms ($R_{\rm H-H}$) in a.u.

The CC and DUCC results for various geometries of the H8 system are collected
in Table 2. In this case, we also used active spaces defined by the two
highest occupied orbitals (orbitals 3 and 4) and the two lowest virtual
orbitals (orbitals 5 and 6). For the equilibrium geometry, including the
triple and quadruple excitations in the external cluster operator, results in
the near-FCI quality of the DUCC results. The accuracy of the DUCC-SD
formalism is also satisfactory (less than 0.7 miliHartree error with respect
to the FCI result), showing the effectiveness of DUCC formalism in compressing
the dynamical correlation effects even in the case when a simple form of
external amplitudes is invoked. In analogy to the H6 case, one can observe the
variational breakdown of the CCSD and CCSDT results for the larger H-H
separations and restoration of the variational character by the corresponding
DUCC formulations.

The DUCC case of downfolding is active-space specific (the so-called cluster
amplitudes universal problem discussed in Ref.Kowalski (2021)). However, to
explore to what extent DUCC downfolding formalism depends on the active space
definition, in Table 3, we collected DUCC-CCSDTQ energies for the H8 model
with $R_{\rm H-H}$=2.0 a.u. and three model spaces defined by $\\{3,4,5,6\\}$,
$\\{2,3,6,7\\}$, and $\\{1,2,7,8\\}$ orbitals. One can see that although the
DUCC-CCSDTQ energies are not equal, the energy discrepancies between results
corresponding to various active spaces do not exceed 0.016 milliHartree.

FCI | DUCC-CCSDTQ | DUCC-CCSDTQ | DUCC-CCSDTQ
---|---|---|---
| ($R=\\{3,4\\}$,$S=\\{5,6\\}$) | ($R=\\{2,3\\}$,$S=\\{6,7\\}$) | ($R=\\{1,2\\}$,$S=\\{7,8\\}$)
-4.286011 | -4.286005 | -4.285865 | -4.285853
Table 3: The DUCC-CCSDTQ results for the H8 model ($R_{\rm H-H}$=2.0 a.u.)
were obtained with the STO-3G basis set for various choices of active spaces.

The effect of the size of the active space is discussed in Table 4 where we
collated DUCC-CCSDTQ results obtained for various geometries of H8 and large
$\\{2,3,4,5,6,7\\}$-generated active space. As expected, the increase in the
active space size results in more accurate DUCC-CCSDT energies, especially for
large distances. It is also evident that the DUCC-CCSDTQ results surpass the
CCSDTQ accuracy for $R_{\rm H-H}$=3.0 a.u. due to the fact that higher-than-
quadruple excitations are included in the diagonalization of the downfolded
Hamiltonian in the active space. This behavior reflects the fact discussed in
Ref.Bauman and Kowalski (2022b) that size of the active space can compensate
for the missing correlation effects not included in the external cluster
operators.

$R_{\rm H-H}$ | FCI | CCSDTQ | DUCC-CCSDTQ
---|---|---|---
2.00 | -4.286011 | -4.286013 | -4.286008
2.50 | -4.114829 | -4.114795 | -4.114782
3.00 | -3.944748 | -3.943920 | -3.944137
Table 4: The DUCC-CCSDTQ results were obtained for various geometries of the
H8 model in the STO-3G basis set for $\\{2,3,4,5,6,7\\}$-generated active
space (see text for details).

From the point of view of quantum computing application, it is essential to
evaluate the performance of techniques to approximate the many-body forms of
the downfolded Hamiltonians and the potential accuracies of quantum solvers.
To this end, we will analyze the finite rank expansion for Eq.(7) and variant
of the Connected Methods Expansion (CMX) Horn and Weinstein (1984); Cioslowski
(1987); Knowles (1987); Mancini, Zhou, and Meier (1994); Noga, Szabados, and
Surján (2002) based on the Peeters-Devreese-Soldatov (PDS) functional Peeters
and Devreese (1984); Soldatov (1995), which has been recently explored in the
context of quantum computing.Kowalski and Peng (2020); Peng and Kowalski
(2021); Claudino  _et al._ (2021) In this paper we will apply low-rank PDS
formulations PDS(3) and PDS(4) (for details we refer the reader to
Refs.Peeters and Devreese (1984); Soldatov (1995)) to identify ground-state
energies of the downfolded Hamiltonians. It is worth noting that the low ranks
of the PDS approach can be effectively implemented on quantum computers. The
PDS results are collated in Table 6, providing Hartree-Fock, complete active
space self-consistent field (using four active electrons distributed over four
active orbitals), and active-space FCI energies. For H6 and H8 models, we used
$\\{2,3,4,5\\}$\- and $\\{4,5,6,7\\}$-generated active spaces, respectively.
Before discussing the PDS results, we should stress the efficiency of the
downfolding procedures (here illustrated in the example of the DUCC-CCSDTQ
approach) in capturing the out-of-active-space correlation effect. This is
best illustrated by comparing DUCC-CCSDTQ vs. CASSCF(4,4) and Active-space FCI
energies. Despite using the same active space definitions, the CASSCF(4,4) and
Active-space FCI energies for all geometries of H6 and H8, in contrast to the
DUCC-CCSDTQ approach, are characterized by significant errors with respect to
the FCI energies. Although, in many cases, quantum simulations are performed
for small dimensionality active spaces using bare Hamiltonians, the quality of
the results can be significantly improved without a significant increase in
quantum computing resources by using a downfolded form of the Hamiltonian. As
seen from Table 5, the PDS(3) can provide much better quality results than
active-space FCI. The PDS(4) can further refine the accuracies of the PDS(3)
approach reducing the errors with respect to exact DUCC-CCSDTQ energies to
within 0.7 milliHartree for H6 and H8 model systems.

Method | H6 | H6 | H8 | H8
---|---|---|---|---
| ($R_{\rm H-H}=2.0$ a.u) | ($R_{\rm H-H}=3.0$ a.u) | ($R_{\rm H-H}=2.0$ a.u) | ($R_{\rm H-H}=3.0$ a.u)
HF | -3.105850 | -2.675432 | -4.138199 | -3.572347
CASSCF(4,4) | -3.175370 | -2.856832 | -4.205528 | -3.699677
Active-space FCI | -3.166938 | -2.802092 | -4.190602 | -3.665605
FCI | -3.217699 | -2.957646 | -4.286011 | -3.944748
DUCC-CCSDTQ | -3.217697 | -2.957384 | -4.286005 | -3.943614
PDS(3) | -3.214888 | -2.953067 | -4.283332 | -3.941223
PDS(4) | -3.217234 | -2.956712 | -4.285622 | -3.943349
Table 5: Comparison of the CASSCF(4,4), active-space FCI, DUCC-CCSDTQ, and
PDS(3)/PDS(4) energies for H6 and H8 model systems. The PDS(3)/PDS(4)
approaches were applied to evaluate the ground-state energy of the DUCC-CCSDTQ
effective Hamiltonians. The $\\{2,3,4,5\\}$\- and $\\{4,5,6,7\\}$-generated
active spaces were used for H6 and H8 systems, respectively.

In the last part of our discussion, we analyze the accuracies of the finite
rank (${\rm Max}_{R}$) approximations for downfolded Hamiltonians. The results
for H6 and H8 models are collected in Table 6. The convergence of the
commutator expansions is illustrated in the example of the 10th-rank
commutator expansion. In all cases discussed in Table 6 ${\rm Max}_{R}$=10
approximation reproduces virtually the exact DUCC-CCSDTQ energies. In
practical applications based on the many-body form of the downfolded
Hamiltonian, only low-rank commutator expansions (${\rm Max}_{R}$=1,…,4) are
numerically feasible (the ${\rm Max}_{R}$=0 corresponds to the active-space
FCI results). One can observe that for weakly correlated situations ($R_{\rm
H-H}$=2.0 a.u.) ${\rm Max}_{R}$=3 provides a satisfactory approximation of the
exact DUCC-CCSDTQ energies. For strongly correlated case ($R_{\rm H-H}$=3.0
a.u.), the inclusion of the 4-th rank commutators (${\rm Max}_{R}$=4) is
needed. In practical application, however, all expansions are based on the
mixing of ${\rm Max}_{R}=n$ contributions with $n+1$-rank commutators stemming
from the Fock matrix terms to reinstate the so-called perturbative balance
(see Refs. for the discussion). For example, ${\rm Max}_{R}$=1 case epitomizes
a situation where the perturbative balance is violated, and non-variational
energies can be obtained. Therefore for the strongly correlated cases, we
recommend the expansions based on the inclusion of ${\rm Max}_{R}$=3 or ${\rm
Max}_{R}$=4 terms with the Fock-operator-dependent terms originating in the
4-th and 5-th commutators, respectively.

${\rm Max_{R}}$/Method | H6 | H6 | H8 | H8
---|---|---|---|---
| ($R_{\rm H-H}=2.0$ a.u) | ($R_{\rm H-H}=3.0$ a.u) | ($R_{\rm H-H}=2.0$ a.u) | ($R_{\rm H-H}=3.0$ a.u)
${\rm Max_{R}}$=0 | -3.166938 | -2.802092 | -4.190602 | -3.665605
${\rm Max_{R}}$=1 | -3.269110 | -3.116145 | -4.382423 | -4.228605
${\rm Max_{R}}$=2 | -3.218732 | -2.976207 | -4.288761 | -3.986313
${\rm Max_{R}}$=3 | -3.217344 | -2.949796 | -4.285044 | -3.927241
${\rm Max_{R}}$=4 | -3.217693 | -2.956814 | -4.285985 | -3.941744
${\rm Max_{R}}$=5 | -3.217699 | -2.957632 | -4.286012 | -3.944383
${\rm Max_{R}}$=10 | 3.217697 | -2.957384 | -4.286005 | -3.943615
DUCC-CCSDTQ | -3.217697 | -2.957384 | -4.286005 | -3.943614
FCI | -3.217699 | -2.957646 | -4.286011 | -3.944748
Table 6: Comparison of energies obtained with finite Commutator Expansion
(CCSDTQ based) for the downfolded Hamiltonians and Exact Downfolding(DUCC-
CCSDTQ) for H6 and H8 models in the STO-3G basis set.

##  V Conclusion

A series of calculations were conducted to examine the impact of
approximations made on external cluster amplitudes on CC downfolded energies.
Simple model systems, H6 and H8 linear chains, were utilized to continuously
vary the extent of correlation effects from weakly to strongly correlated
regimes (with $R_{\rm H-H}$ ranging from 2.0 to 3.0 a.u.). The results showed
that while the external cluster amplitudes from SR-CCSD calculations were
satisfactory for the weakly correlated situation, for the strongly correlated
case, the effect of triply and quadruply excited external clusters could no
longer be neglected. The downfolding procedure acted as stabilizers and could
restore the variational character of energies despite the fact that external
amplitudes were obtained from SR-CC calculations that suffer from the
variational collapse. Furthermore, the downfolded energies obtained for
various active spaces on the H8 systems (including those that did not include
essential correlation effects) had only small discrepancies, demonstrating the
approximate invariance of downfolding energies for various active spaces. This
was the case for all SES-CC types active spaces for commutative SR-CC
formulations. Additionally, it was shown that the downfolded Hamiltonians
could be effectively diagonalized with low-order PDS formulations for both
weakly and strongly correlated regimes.

The assessment of the impact of the maximum rank of the commutator expansion
on the precision of downfolded energies is an integral aspect of the analysis
associated with practical applications of CC downfolding procedures. Our
findings demonstrate that an increase in the degree of correlation effects
necessitates the inclusion of higher rank commutators. In particular, for
H6/H8 $R_{\rm H-H}=2.0$ a.u., the inclusion of single, double, and triple
commutators is sufficient to achieve a favorable agreement with the FCI
energies. However, for $R_{\rm H-H}=3.0$ a.u., higher-order commutators
(quadruple/pentuple) must be incorporated into the approximation. It should be
noted, however, that for the $R_{\rm H-H}=3.0$ a.u. scenario, an active space
that can distinguish between static and dynamical correlation effects cannot
be constructed (i.e. all orbitals must be considered active). Nonetheless, as
Table (6) demonstrates, the corresponding commutator-rank expansion rapidly
converges to the FCI energies.

##  VI Acknowledgement

This work was supported by the Quantum Science Center (QSC), a National
Quantum Information Science Research Center of the U.S. Department of Energy
(under FWP 76213) and by “Embedding QC into Many-body Frameworks for Strongly
Correlated Molecular and Materials Systems” project, which is funded by the
U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences,
the Division of Chemical Sciences, Geosciences, and Biosciences (under FWP
72689). This work used resources from the Pacific Northwest National
Laboratory (PNNL). PNNL is operated by Battelle for the U.S. Department of
Energy under Contract DE-AC05-76RL01830.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts of interest to declare.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.

## References

  * Takeshita  _et al._ (2020) T. Takeshita, N. C. Rubin, Z. Jiang, E. Lee, R. Babbush, and J. R. McClean, Phys. Rev. X 10, 011004 (2020).
  * Şahinoğlu and Somma (2021) B. Şahinoğlu and R. D. Somma, npj Quantum Information 7, 119 (2021).
  * Liu, Li, and Yang (2022) J. Liu, Z. Li, and J. Yang, Journal of Chemical Theory and Computation 18, 4795 (2022).
  * Huang, Li, and Evangelista (2023) R. Huang, C. Li, and F. A. Evangelista, PRX Quantum 4, 020313 (2023).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, Nat. Commun. 5, 4213 (2014).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016).
  * Romero  _et al._ (2018) J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Quantum Sci. Technol. 4, 014008 (2018).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * Izmaylov  _et al._ (2019) A. F. Izmaylov, T.-C. Yen, R. A. Lang, and V. Verteletskyi, J. Chem. Theory Comput. 16, 190 (2019).
  * Lang, Ryabinkin, and Izmaylov (2021) R. A. Lang, I. G. Ryabinkin, and A. F. Izmaylov, JJ. Chem. Theory Comput. 17, 66 (2021).
  * Grimsley  _et al._ (2019a) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nat. Commun. 10, 1 (2019a).
  * Grimsley  _et al._ (2019b) H. R. Grimsley, D. Claudino, S. E. Economou, E. Barnes, and N. J. Mayhall, J. Chem. Theory Comput. 16, 1 (2019b).
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Reviews of Modern Physics 92, 015003 (2020).
  * Kirby and Love (2021) W. M. Kirby and P. J. Love, Phys. Rev. Lett. 127, 110503 (2021).
  * Tilly  _et al._ (2022) J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, _et al._ , Physics Reports 986, 1 (2022).
  * Coester (1958) F. Coester, Nucl. Phys. 7, 421 (1958).
  * Coester and Kummel (1960) F. Coester and H. Kummel, Nucl. Phys. 17, 477 (1960).
  * Čížek (1966) J. Čížek, J. Chem. Phys. 45, 4256 (1966).
  * Paldus, Čížek, and Shavitt (1972) J. Paldus, J. Čížek, and I. Shavitt, Phys. Rev. A 5, 50 (1972).
  * Mukherjee, Moitra, and Mukhopadhyay (1975) D. Mukherjee, R. K. Moitra, and A. Mukhopadhyay, Mol. Phys. 30, 1861 (1975).
  * Adams, Jankowski, and Paldus (1979) B. Adams, K. Jankowski, and J. Paldus, Chem. Phys. Lett. 67, 144 (1979).
  * Purvis and Bartlett (1982) G. D. Purvis and R. J. Bartlett, J. Chem. Phys. 76, 1910 (1982).
  * Jeziorski and Monkhorst (1981) B. Jeziorski and H. J. Monkhorst, Phys. Rev. A 24, 1668 (1981).
  * Lyakh  _et al._ (2012) D. I. Lyakh, M. Musiał, V. F. Lotrich, and R. J. Bartlett, Chem. Rev. 112, 182 (2012).
  * Bauman  _et al._ (2019) N. P. Bauman, E. J. Bylaska, S. Krishnamoorthy, G. H. Low, N. Wiebe, C. E. Granade, M. Roetteler, M. Troyer, and K. Kowalski, J. Chem. Phys. 151, 014107 (2019).
  * Kowalski (2018) K. Kowalski, J. Chem. Phys. 148, 094104 (2018).
  * Kowalski (2021) K. Kowalski, Phys. Rev. A 104, 032804 (2021).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Scuseria and Schaefer (1988) G. E. Scuseria and H. F. Schaefer, Chem. Phys. Lett. 152, 382 (1988).
  * Noga and Bartlett (1987) J. Noga and R. J. Bartlett, J. Chem. Phys. 86, 7041 (1987).
  * Noga and Bartlett (1988) J. Noga and R. J. Bartlett, J. Chem. Phys. 89, 3401 (1988).
  * Oliphant and Adamowicz (1991a) N. Oliphant and L. Adamowicz, J. Chem. Phys. 95, 6645 (1991a).
  * Piecuch and Adamowicz (1994) P. Piecuch and L. Adamowicz, J. Chem. Phys. 100, 5792 (1994).
  * Kucharski and Bartlett (1991) S. A. Kucharski and R. J. Bartlett, Theor. Chem. Acc. 80, 387 (1991).
  * Kowalski and Bauman (2020) K. Kowalski and N. P. Bauman, J. Chem. Phys. 152, 244127 (2020).
  * Bauman and Kowalski (2022a) N. P. Bauman and K. Kowalski, Materials Theory 6, 1 (2022a).
  * Oliphant and Adamowicz (1991b) N. Oliphant and L. Adamowicz, J. Chem. Phys. 94, 1229 (1991b).
  * Oliphant and Adamowicz (1992) N. Oliphant and L. Adamowicz, J. Chem. Phys. 96, 3739 (1992).
  * Piecuch, Oliphant, and Adamowicz (1993) P. Piecuch, N. Oliphant, and L. Adamowicz, J. Chem. Phys. 99, 1875 (1993).
  * Hehre, Stewart, and Pople (1969) W. J. Hehre, R. F. Stewart, and J. A. Pople, The Journal of Chemical Physics 51, 2657 (1969).
  * Bauman and Kowalski (2022b) N. P. Bauman and K. Kowalski, J. Chem. Phys 156, 094106 (2022b).
  * Horn and Weinstein (1984) D. Horn and M. Weinstein, Phys. Rev. D 30, 1256 (1984).
  * Cioslowski (1987) J. Cioslowski, Phys. Rev. Lett. 58, 83 (1987).
  * Knowles (1987) P. J. Knowles, Chem. Phys. Lett. 134, 512 (1987).
  * Mancini, Zhou, and Meier (1994) J. D. Mancini, Y. Zhou, and P. F. Meier, International Journal of Quantum Chemistry 50, 101 (1994).
  * Noga, Szabados, and Surján (2002) J. Noga, A. Szabados, and P. R. Surján, International Journal of Molecular Sciences 3, 508 (2002).
  * Peeters and Devreese (1984) F. Peeters and J. Devreese, Journal of Physics A: Mathematical and General 17, 625 (1984).
  * Soldatov (1995) A. Soldatov, International Journal of Modern Physics B 9, 2899 (1995).
  * Kowalski and Peng (2020) K. Kowalski and B. Peng, The Journal of Chemical Physics 153, 201102 (2020).
  * Peng and Kowalski (2021) B. Peng and K. Kowalski, Quantum 5, 473 (2021).
  * Claudino  _et al._ (2021) D. Claudino, B. Peng, N. P. Bauman, K. Kowalski, and T. S. Humble, Quantum Science and Technology 6, 034012 (2021).