# A Fast Numerical Solver of Quantum-inspired
Ising Optimization Problems

Langyu Li, and Yu Pan L. Li and Y. Pan are with the College of Control Science
and Engineering, Zhejiang University, Hangzhou 310027, China.
Email address<EMAIL_ADDRESS>
###### Abstract

Quantum annealers, coherent Ising machines and digital Ising machines for
solving quantum-inspired optimization problems have been developing rapidly
due to their near-term applications. The numerical solvers of the digital
Ising machines are based on traditional computing devices. In this work, we
propose a fast and efficient solver for the Ising optimization problems. The
algorithm consists of a pruning method that exploits the graph information of
the Ising model to reduce the computational complexity, and a domain selection
method which introduces significant acceleration by relaxing the discrete
feasible domain into a continuous one to incorporate the efficient gradient
descent method. The experiment results show that our solver can be an order of
magnitude faster than the classical solver, and at least two times faster than
the quantum-inspired annealers including the simulated quantum annealing on
the benchmark problems. With more relaxed requirements on hardware and lower
cost than quantum annealing, the proposed solver has the potential for near-
term application in solving challenging optimization problems as well as
serving as a benchmark for evaluating the advantage of quantum devices.

###### Index Terms:

Ising optimization, Ising machine, quadratic unconstrained binary
optimization, quantum annealing.

##  I Introduction

As a fundamental concept in the realm of quantum computing and optimization,
the Ising problem has garnered significant attention in recent years. The
Ising problem, originally conceived in statistical mechanics [1], has found
applications in quantum computing [2] which is mainly due to its profound
potential in solving large-scale combinatorial optimization problems and
enhance the performance of machine learning [3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15].

Figure 1: Illustration of our solver. (a) Our solver first applies pruning to
simplify the optimization problem by utilizing the graph information of the
Ising model. (b) Then the domain selection method is applied to find the
approximate optimal solutions. As illustrated in (b1)-(b3), the original
problem is firstly transformed into an $L_{1}$-norm maximization problem with
discrete feasible domains. Then we relax the feasible domain to a continuous
one and use the gradient descent optimizer to find the approximate optimal
domain. (c) Finally, the solution of the primal problem is retrieved by
combining the results of domain selection and pruning.

Both quantum and classical solvers can be used to solve the Ising optimization
problems. Quantum solvers mainly include quantum annealers (QAs), and coherent
Ising machines (CIMs) based on superconducting and optical devices,
respectively. Grounded on quantum physics and quantum adiabatic theorem,
quantum solvers leverage quantum properties such as superposition and
entanglement to explore the energy landscapes of complex systems [16].
Pioneered by D-Wave, QAs have garnered significant attention due to their
potential to outperform classical computing in solving optimization tasks [17,
18]. By employing degenerate optical parametric oscillators, CIMs exploit
optical interference phenomena to encode and solve complex Ising problems [19,
20, 21]. Although QAs offer unique capabilities with the inherent quantum
parallelism, their severe operating constraints such as low temperature and
size make it difficult to deploy on small intelligent devices [22, 23]. CIMs
face similar limitations [24].In contrast, quantum-inspired classical solvers
are much cheaper and smaller. Despite the lack of quantum advantages inherent
in quantum devices, current classical solvers can achieve speed and solution
quality close to that of quantum solvers. Therefore, digital solvers based on
traditional devices are more likely to be implemented for realistic
applications in the near term. Besides, classical solvers are necessary
benchmarks for understanding and measuring the advantage of quantum solvers.

Classical solvers mainly refer to digital Ising machines (DIMs). Although
rooted in classical computing, DIMs have shown potential in solving complex
Ising optimization problem [25, 26, 27, 28]. DIMs employ digital hardware to
implement the optimization algorithms. The architectures of digital Ising
machines can be categorized into non-von Neumann architecture and von Neumann.
The former often requires an application-specific integrated circuit (ASIC)
that uses a large number of registers to simulate the spins in the Ising
model, and the update of the register states is controlled by the simulated
annealing algorithm [29, 30, 31]. The latter uses the general-purpose
computing devices such as CPU, GPU, and field-programmable gate arrays (FPGAs)
without physically simulating the spins of the Ising model. In particular,
recent works based on non-Von Neumann architectures [31, 26, 28] have not been
able to shown a significant advantage over the von Neumann architecture [32].

Among the methods based on von Neumann architecture, the operator-based
methods [33, 34] tend to have high computational complexity and are not
suitable for large-scale problems. Meanwhile, the accuracy and speed of the
other numerical solvers [35, 36] still have room for improvement. In this
work, we propose a fast and scalable solver which demonstrates an advantage in
speed and accuracy. This solver first applies pruning to simplify and reduce
the dimensions of the Ising optimization problem by pre-determining the
directions of certain spins. After this preprocessing, the domain selection
method will solve the simplified problem to find the approximate optimal
solutions. In contrast to the previous works, the main innovation is the
fusion of the discrete programming method with the gradient descent. To be
more specific, we transform the Ising problem into an $L_{1}$-norm
maximization problem, which is then decomposed into several sub-linear
programming problems. The feasible sets between these linear programming
problems are not continuous, and the optimal solution of each sub-problem is
easy to find. Therefore, we find a connection between the energy of the
feasible solution and the feasible domain, and design loss functions weighted
by the energy of the feasible solution. By making the feasible domain
continuous, gradient descent can be used to quickly find a feasible domain
containing a near-optimal solution, which is the optimal solution in this
domain. In the last, retrieve the solutions of the primal Ising problem from
the solution of domain selection with the information recorded in pruning.

A series of experiments have been conducted to evaluate the performance of our
solver. In particular, our solver has been extensively compared with Gurobi,
D-Wave classical steepest descent solver, simulated annealing and simulated
quantum annealing. Our solver is at least $10$ times faster than Gurobi, which
is the state-of-the-art solver based on operation programming. Our solver is
also significantly faster than other classical numerical methods and quantum
simulated annealing. The comparisons are made on randomly generated samples,
as well as some of the most famous NP-hard optimization problems such as
maximum cut. The computational complexity has also been verified with the
experiment data, which demonstrates the outstanding scalability of our solver
as the problem size scales up.

This paper is organized as follows. Sec. II introduces the mathematical model
of Ising optimization problem. Sec. III describes our solver, with Sec. III-A
and Sec. III-B introducing the pruning and domain selection algorithms,
respectively. Sec. IV summarizes the experiment results and comparison. Sec. V
concludes this work.

##  II Ising problem

The Ising model was first formulated by Ernst Ising in 1925 [1], which was
initially conceived as a simple model for describing magnetic interactions in
crystalline systems explaining the behavior of spins interacting with their
neighbors. The Hamiltonian of the Ising model is defined by

$\displaystyle
H=\sum_{i>j}\widehat{J}_{ij}\sigma_{i}^{z}\sigma_{j}^{z}+\sum_{i}\widehat{h}_{i}\sigma_{i}^{z},$
(1)

where $\sigma_{i}^{z}$ is the Pauli-$Z$ operator acting on the $i$-th qubit.
The coupling amplitude $\widehat{J}_{ij}\in\mathbb{R}$ is a scalar to describe
the strength of the coupling between the qubits $i$ and $j$, and the
coefficient $\widehat{h}_{i}\in\mathbb{R}$ describes the external field
strength applied on the $i$-th qubit. The qubits are also called spins since
each qubit has two states $|0\rangle$ and $|1\rangle$ corresponding to spin up
and spin down, respectively. Fig. 2 shows a lattice Ising model.

Figure 2: A lattice Ising model without external field.

The state of the Ising model can be described by a state vector
$|\psi\rangle\in\\{|0\rangle,|1\rangle\\}^{\otimes d}$. For a $d$-qubit Ising
model with the state vector $|\psi\rangle$, its energy, or called landscape,
is given as

$\displaystyle\langle
H\rangle=\langle\psi|H|\psi\rangle,\quad|\psi\rangle\in\\{|0\rangle,|1\rangle\\}^{\otimes
d}.$ (2)

Eq. (2) has an equivalent scalar form as

$\displaystyle\langle H\rangle$ $\displaystyle\cong E(\bm{s}),$ (3)
$\displaystyle E(\bm{s})$
$\displaystyle=\bm{s}^{\top}\bm{J}\bm{s}+\bm{h}^{\top}\bm{s},\quad\bm{s}\in\\{-1,1\\}^{d},$
(4)

where $\bm{s}$ is a $d$-dimensional vector called spin configuration with
respect to the state vector $|\psi\rangle$. The $i$-th entry of $\bm{s}$ is
$+1$ if the $i$-th spin is at the state $|0\rangle$, and $-1$ if the spin
state is $|1\rangle$. Coupling matrix $\bm{J}\in\mathbb{R}^{d\times d}$ is
defined by $J_{ij}=J_{ji}=\frac{1}{2}\widehat{J}_{ij}$, which is real
symmetric with $J_{ii}=0$. Field vector $\bm{h}\in\mathbb{R}^{d}$ is defined
by $h_{i}=\widehat{h}_{i}$. The flip operator on the $i$-th spin is defined as

$\displaystyle
X_{i}=\operatorname{diag}(1,\cdots,1,\underset{i\text{-th}}{-1},1,\cdots,1).$
(5)

We use the notation $X_{i,j}=X_{i}X_{j}$ since $X_{i}$ and $X_{j}$ commute.

The Ising optimization problem refers to the task of finding the optimal
configuration of spins within a system to minimize or maximize the given
objective function. In particular, most of the Ising optimization problem can
be formulated as finding the ground-state energy of the Ising model.
Mathematically, we can define this primal problem as follows

$\displaystyle\min_{\bm{s}}\quad$ $\displaystyle E(\bm{s}),$ (6) s.t.
$\displaystyle\bm{s}\in\\{-1,1\\}^{d}.$ (7)

$\bm{s}^{*}$ denotes the ground state vector or optimal state vector, which is
written as

$\displaystyle\bm{s}^{*}=\arg\min_{\bm{s}}\quad$ $\displaystyle
E(\bm{s}),\quad\bm{s}\in\\{-1,1\\}^{d}.$ (8)

##  III Pruning and Domain Selection

In this section, we will introduce the numerical solver which consists of two
stages, namely pruning and domain selection. Pruning can be considered as a
preprocessing step of the primal Ising problem based on the graph information,
and domain selection is an efficient optimization method to produce the near-
optimal solution.

###  III-A Pruning

Pruning is the process of pre-determining the optimal directions of certain
spins. For an arbitrary state $\bm{s}$, its energy difference between an
arbitrarily flipped spin configuration $\bm{s}^{\prime}$ is calculated as

$\displaystyle
E(\bm{s}^{\prime})-E(\bm{s})=E(X_{\mathcal{F}}\bm{s})-E(\bm{s})$ (9)
$\displaystyle=$
$\displaystyle\bm{s}^{\top}(X_{\mathcal{F}}\bm{J}X_{\mathcal{F}}-\bm{J})\bm{s}+(\bm{h}^{\top}X_{\mathcal{F}}-\bm{h}^{\top})\bm{s}$
(10) $\displaystyle=$
$\displaystyle-4\sum_{i\in\mathcal{F}}s_{i}\sum_{j=1}^{d}J_{ij}s_{j}-2\sum_{i\in\mathcal{F}}h_{i}s_{i},$
(11)

where $X_{\mathcal{F}}=\prod_{i\in\mathcal{F}}X_{i}$. For the $i$-th spin, if
there holds $2\sum_{j}|J_{ji}|<|h_{i}|$, we must have

$\displaystyle s^{*}_{i}=-\operatorname{sign}(h_{i}),$ (12)

since the energy of the configuration will definitely increase if $s^{*}_{i}$
is flipped

$\displaystyle E(X_{i}\bm{s})-E(\bm{s})$ (13) $\displaystyle=$
$\displaystyle-4s_{i}\sum_{j=1}^{d}J_{ij}s_{j}-2h_{i}s_{i}$ (14)
$\displaystyle=$
$\displaystyle-4s_{i}\sum_{j=1}^{d}J_{ij}s_{j}+2\operatorname{sign}(h_{i})h_{i}$
(15) $\displaystyle\geq$ $\displaystyle
2|h_{i}|-4\big{|}s_{i}\sum_{j=1}^{d}J_{ij}s_{j}\big{|}>0.$ (16)

This class of spins are named as determined spins. Similarly, if the $i$-th
spin is coupled to the $j$-th spin only and there holds
$2|J_{ij}|\geq|h_{i}|$, we must have

$\displaystyle s^{*}_{i}=-s^{*}_{j}\operatorname{sign}(J_{ij}).$ (17)

This class of spins are named as semi-determined spins since their directions
are solely dependent on the neighboring spins.

Therefore, given an Ising problem, the spins can be divided into determined
spins $\mathcal{D}$, semi-determined spins $\mathcal{S}$, and other spins
$\mathcal{G}$. Without loss of generality, we can write the simplified
optimization problem after a single pruning as

$\displaystyle\min_{\bm{s}_{\mathcal{G}}}\quad$
$\displaystyle\bm{s}^{\top}_{\mathcal{G}}\bm{J}_{\mathcal{G}}\bm{s}_{\mathcal{G}}+\bm{h}^{\top}_{\mathcal{G}}\bm{s}_{\mathcal{G}}+C_{\mathcal{D},\mathcal{S}},$
(18) s.t. $\displaystyle\bm{s}_{\mathcal{G}}\in\\{-1,1\\}^{|\mathcal{G}|},$
(19)

where

$\displaystyle(\bm{J}_{\mathcal{G}})_{ij}=$ $\displaystyle
J_{ij},\quad\text{for $i,j\in\mathcal{G}$},$ (20)
$\displaystyle(\bm{h}_{\mathcal{G}})_{i}=$ $\displaystyle
h_{i}-2\sum_{j\in\mathcal{D}}J_{ij}\operatorname{sign}(h_{j})-\sum_{k\in\mathcal{S}}\operatorname{sign}(J_{ik})h_{k}$
(21)
$\displaystyle+2\sum_{k\in\mathcal{S}}\sum_{j\in\mathcal{D}}J_{kj}\operatorname{sign}(J_{ik}h_{j}),\quad\text{for
$i\in\mathcal{G}$},$ (22) $\displaystyle C_{\mathcal{D},\mathcal{S}}=$
$\displaystyle-\sum_{i\in\mathcal{D}}|h_{i}|-2\sum_{j\in\mathcal{G}}\sum_{k\in\mathcal{S}}|J_{jk}|+\sum_{i,j\in\mathcal{D}}J_{ij}\operatorname{sign}(h_{i}h_{j}).$
(23)

The pruning algorithm can then be defined recursively according to Alg. 1. An
example of pruning is shown in Fig. 3.

input : $d$-spin Ising problem with coupling matrix $\bm{J}$ and field vector
$\bm{h}$

output : A $d^{\prime}$-spin Ising problem with matrix $\bm{J}^{\prime}$ and
vector $\bm{h}^{\prime}$,$d^{\prime}\leq d$

1

2while  _$\mathcal{S}\cup\mathcal{D}\neq\emptyset$_ do

3

4 Select the determined spins into $\mathcal{D}$;

5 Select the semi-determined spins into $\mathcal{S}$;

6 Divide the current spins into $\mathcal{G},\mathcal{S},\mathcal{D}$;

7 Redefine the current Ising problem by Eq. (18) $\sim$ (23).

8 end while

return The simplified Ising problem with $\bm{J}^{\prime}$ and
$\bm{h}^{\prime}$.

Algorithm 1 Pruning Figure 3: A $6$-spin example of the pruning algorithm.
According to Eq. (12), the spin direction of $s_{4}$ is down. After
determining $s_{4}=-1$, the primal optimization problem is divided into two
sub-problems with updated parameters. Then one more spin direction $s_{6}=1$
can be determined in the bottom triangle.

###  III-B Domain Selection

The domain selection method first converts the primal problem into an
$L_{1}$-norm maximization problem whose feasible domain is not continuous.
Then the feasible domain is relaxed to a continuous one to facilitate
gradient-based optimization method.

By introducing another variable we can transform the primal problem into

$\displaystyle\min_{\bm{s}_{x},\bm{s}_{y}}$
$\displaystyle\quad\bm{s}_{x}^{\top}\bm{J}\bm{s}_{y}+\bm{h}^{\top}\bm{s}_{y},$
(24) s.t. $\displaystyle\quad\bm{s}_{x}=\bm{s}_{y},$ (25)
$\displaystyle\quad\bm{s}_{x},\bm{s}_{y}\in\\{-1,1\\}^{d}.$ (26)

The following metric function $D(\bm{s}_{x},\bm{s}_{y})$

$\displaystyle D(\bm{s}_{x},\bm{s}_{y})=d-\bm{s}_{x}^{\top}\bm{s}_{y}$ (27)

is used to measure the distance between $\bm{s}_{x}$ and $\bm{s}_{y}$. Thus
the augmented Lagrange formulation of this two-variable problem is given by

$\displaystyle\min_{\bm{s}_{x},\bm{s}_{y}}$
$\displaystyle\quad\bm{s}_{x}^{\top}\bm{J}\bm{s}_{y}+\bm{h}^{\top}\bm{s}_{y}+\mu(d-\bm{s}_{x}^{\top}\bm{s}_{y}),$
(28) s.t. $\displaystyle\quad\bm{s}_{x},\bm{s}_{y}\in\\{-1,1\\}^{d},$ (29)

where $\mu>0$ is a weight scalar. The optimal solution to this problem can be
calculated as

$\displaystyle\min_{\bm{s}_{x},\bm{s}_{y}}\quad\bm{s}_{x}^{\top}\bm{J}\bm{s}_{y}+\bm{h}^{\top}\bm{s}_{y}-\mu\bm{s}_{x}^{\top}\bm{s}_{y}$
(30) $\displaystyle=$
$\displaystyle\min_{\bm{s}_{x},\bm{s}_{y}}\quad(\bm{s}_{x}^{\top}\bm{J}+\bm{h}^{\top}-\mu\bm{s}_{x}^{\top})\bm{s}_{y}$
(31) $\displaystyle=$
$\displaystyle\min_{\bm{s}_{x}}\min_{\bm{s}_{y}}\quad(\bm{s}_{x}^{\top}\bm{J}+\bm{h}^{\top}-\mu\bm{s}_{x}^{\top})\bm{s}_{y}$
(32) $\displaystyle=$
$\displaystyle\min_{\bm{s}_{x}}\quad-\|\bm{J}\bm{s}_{x}+\bm{h}-\mu\bm{s}_{x}\|_{1}$
(33) $\displaystyle=$ $\displaystyle\max_{\bm{s}}\quad\|(\bm{J}-\mu
I)\bm{s}+\bm{h}\|_{1}.$ (34)

Note that in Eq (32)$\sim$(33) the optimal $\bm{s}_{y}^{*}$ is given by
$(\bm{s}_{y}^{*})_{i}=-\operatorname{sign}((\bm{s}_{x}^{\top}\bm{J}+\bm{h}^{\top}-\mu\bm{s}_{x}^{\top})_{i})$
and $\|\cdot\|_{1}$ is the $L_{1}$ norm. A large enough $\mu$ can guarantee
that the solution to the Augmented Lagrangian problem is the same as the
solution to the primal problem.

The $L_{1}$-norm maximization defined in (34) can then be decomposed into
$2^{d}$ linear programming problems. Then the $L_{1}$-norm optimization
problem has an equivalent form of

$\displaystyle\max_{\bm{w}}$ $\displaystyle\quad\bm{1}^{\top}\bm{w},$ (35)
s.t. $\displaystyle\quad\bm{w}\in\mathcal{C},$ (36)

where $\bm{1}$ is the all-one vector. The feasible set of solutions is the
concatenation of $2^{d}$ convex sets, which is given by

$\displaystyle\mathcal{C}=$
$\displaystyle\bigcup_{i=0}^{2^{d}-1}\mathcal{C}_{i},$ (37)

and each feasible set $\mathcal{C}_{i}$ is constrained by

$\displaystyle\mathcal{C}_{i}=\\{\bm{w}:\bm{w}+X_{\\{i\\}}((\bm{J}-\mu
I)\bm{s}+\bm{h})\leq 0,\bm{s}\in\\{-1,1\\}^{d}\\},$ (38)

with $X_{\\{i\\}}=\prod_{j\in\text{idx}(i)}X_{j}$. $\text{idx}(i)$ returns the
index set of $d$-dimension binary vector whose entry value is equal to 1, and
the binary vector can be seen as a binary number whose value in decimal is
$i$. Consequently, the $L_{1}$-norm optimization problem becomes a feasible
domain search for the $\mathcal{C}_{i}$ where the optimal solution locates.
Finding such an optimal value in a certain feasible domain is a mixed-integer
linear programming problem, which has exponential complexity in principle
[37]. In our case, this problem is relatively easy to solve since the local
optimal solution of $\mathcal{C}_{i}$ can be obtained as

$\displaystyle\widetilde{\bm{s}}^{*}=-\operatorname{sign}(X_{\\{i\\}}(\bm{J}-\mu
I)\bm{1})=\operatorname{sign}(X_{\\{i\\}}\bm{1}).$ (39)

Therefore, all we need to do is to find the optimal domain
$\mathcal{C}_{i}^{*}$ if we encode the local optimal solutions into the
feasible domains such that different feasible domains have different costs.
The discrete $X_{\\{i\\}}$ can be relaxed by multiplying a factor of matrix
function $\sin(\bm{\theta})$ as

$\displaystyle\max_{\bm{w}}$ $\displaystyle\quad\bm{1}^{\top}\bm{w},$ (40)
s.t. $\displaystyle\quad\bm{w}+\sin{(\bm{\theta})}((\bm{J}-\mu
I)\widetilde{\bm{s}}^{*}+\bm{h})\leq 0.$ (41)

where
$\bm{\theta}=\operatorname{diag}\\{\theta_{1},\theta_{2},\cdots,\theta_{d}\\}$.
Such a relaxation transforms the original discrete parameter space into a
continuous parameter space. By incorporating the constraints into the
optimization target, the loss function for domain selection finally becomes

$\displaystyle\mathcal{L}(\bm{\theta})=\bm{1}^{\top}\sin{(\bm{\theta})}((\bm{J}-\mu
I)\widetilde{\bm{s}}^{*}+\bm{h}).$ (42)

With the loss function defined in Eq. (42), the domain selection can be
efficiently solved with any gradient descent method by optimizing
$\bm{\theta}$.

##  IV Numerical Results

###  IV-A Pruning Rate

To evaluate the performance of the pruning algorithm, the pruning rate $\eta$
is defined as

$\displaystyle\eta:=1-\frac{d_{\mathcal{G}}}{d},$ (43)

where $d_{\mathcal{G}}=|\mathcal{G}|$ is the amount of spins that are not
pruned. Let $\mathcal{U},\mathcal{N},\mathcal{B}$ denote uniform, Gaussian,
and uniformly binary distributions, respectively. For a $d$-qubit $k$-regular
Ising optimization problem with $k\leq
d-1,\widehat{J}_{ij}\sim\mathcal{U}(-\alpha,\alpha)$ and
$\widehat{h}_{i}\sim\mathcal{U}(-\beta,\beta)$, the pruning rate $\eta_{1}$ of
the first round can be estimated as

$\displaystyle\eta_{1}$
$\displaystyle=\Pr(|\widehat{h}_{i}|>\sum_{j}|\widehat{J}_{ij}|)$ (44)
$\displaystyle\approx\Pr\big{(}b>a:a\sim\mathcal{N}\big{(}\frac{k\alpha}{2},\frac{k\alpha^{2}}{12}\big{)},b\sim\mathcal{U}(0,\beta)\big{)}$
(45)
$\displaystyle=\int_{0}^{\beta}\frac{1}{\beta}~{}\text{d}x\int_{-\infty}^{x}\frac{\sqrt{6}}{\sqrt{\pi
k}\alpha}\exp(-6\frac{(y-\frac{k\alpha}{2})^{2}}{k\alpha^{2}})~{}\text{d}y,$
(46)

where we have used the Gaussian distribution to approximate the Irwin–Hall
distribution [38] in Eq. (44) according to the central limit theorem.

###  IV-B Performance

Here we introduce the Hamiltonian ratio $\chi$ as a normalized metric to
evaluate the Hamiltonian error with respect to the ground state [36], which is
defined as

$\displaystyle\chi$
$\displaystyle=\frac{E_{\max}-\widehat{E}}{E_{\max}-E_{\min}},$ (47)

where $\widehat{E}$ is the solution, and $E_{\max}$ and $E_{\min}$ are the
maximal and minimal energies of the Ising model, respectively. Due to the
difficulty in obtaining the exact $E_{\max}$ and $E_{\min}$, $E_{\max}$ and
$E_{\min}$ are estimated by

$\displaystyle\widehat{E}_{\max}=-\widehat{E}_{\min}=\sqrt{d}\|\bm{J}\|_{F}+\|\bm{h}\|_{1},$
(48)

where $\|\cdot\|_{F}$ is the Frobenius norm. The Hamiltonian ratio is thus
approximated by

$\displaystyle\widehat{\chi}$
$\displaystyle=\frac{\widehat{E}_{\max}-\widehat{E}}{\widehat{E}_{\max}-\widehat{E}_{\min}}=\frac{1}{2}-\frac{\widehat{E}}{2(\sqrt{d}\|\bm{J}\|_{F}+\|\bm{h}\|_{1})}.$
(49)

Such an approximation effectively normalizes the Hamiltonian ratio, and
$\widehat{\chi}$ is linearly and positively correlated with $\chi$ as

$\displaystyle\chi=\frac{E_{\max}-(1-2\widehat{\chi})(\sqrt{d}\|\bm{J}\|_{F}+\|\bm{h}\|_{1})}{E_{\max}-E_{\min}}.$
(50)

TABLE I: The performance of pruning and domain selection. PDS and DS refer to
domain selection with and without pruning, respectively. Iterations are the
number of steps taken to reach the convergence. $T$ is the time cost in
seconds.
$\Delta\widehat{\chi}=\widehat{\chi}_{\text{PDS}}-\widehat{\chi}_{\text{DS}}$.

| DS Iterations | $T_{\text{DS}}$ | $\widehat{\chi}_{\text{DS}}$ | PDS Iterations | $T_{\text{PDS}}$ | $\widehat{\chi}_{\text{PDS}}$ | $\eta$ | $\Delta\widehat{\chi}$
---|---|---|---|---|---|---|---|---
Uniform | $89.5\pm 7.4$ | $0.279\pm 0.035$ | $0.9248\pm 0.0037$ | $74.7\pm 12.5$ | $0.056\pm 0.013$ | $0.9261\pm 0.0036$ | $0.949\pm 0.018$ | $0.0014\pm 0.0004$
Gaussian | $94.2\pm 3.4$ | $0.297\pm 0.020$ | $0.9142\pm 0.0052$ | $90.4\pm 5.3$ | $0.246\pm 0.018$ | $0.9145\pm 0.0051$ | $0.290\pm 0.030$ | $0.0003\pm 0.0007$
Binary | $92.7\pm 5.4$ | $0.289\pm 0.028$ | $0.8815\pm 0.0049$ | $93.0\pm 5.2$ | $0.296\pm 0.025$ | $0.8815\pm 0.0049$ | $0.0$ | $0.0000\pm 0.0005$

TABLE II: Details of the experiments in Fig. 4.

| $\mathbb{E}(\widehat{\chi})$ | $\operatorname{Var}(\widehat{\chi})$ | $\mathbb{E}(T)$ | $\operatorname{Var}(T)$ | $\mathbb{E}(\widehat{T}_{\text{TTS}})$ | $\operatorname{Var}(\widehat{T}_{\text{TTS}})$ | $\mathbb{E}(v)$ | $\operatorname{Var}(v)$
---|---|---|---|---|---|---|---|---
Ours | 1.169 | 0.002 | 0.196 | 0.001 | 280.6 | 2183.5 | 6.041 | 0.393
Gurobi | 1.209 | 0.003 | 3.373 | 0.103 | 4641.2 | 377051.4 | 0.361 | 0.001
DWave-SD | 1.149 | 0.002 | 0.319 | 0.002 | 465.9 | 5880.7 | 3.625 | 0.057
OpenJij-SA | 1.115 | 0.002 | 2.127 | 0.039 | 3208.0 | 102109.2 | 0.526 | 0.001
OpenJij-SQA | 0.972 | 0.001 | 2.181 | 0.001 | 3793.2 | 40407.6 | 0.447 | 0.001

TABLE III: Details of the experiments in Fig. 5.

|  |  |  |  | Ours | Gurobi | DWave-SD | OpenJij-SA | OpenJij-SQA
---|---|---|---|---|---|---|---|---|---
Problem | $J_{ij}$ | $h_{i}$ | density |  $k$-regular | $T$ | $\widehat{E}$ | $T$ | $\widehat{E}$ | $T$ | $\widehat{E}$ | $T$ | $\widehat{E}$ | $T$ | $\widehat{E}$
MaxCut-3 | $\\{0,1\\}$ | 0 | 0.003 | 3 | 0.23 | -2312 | 0.16 | -2408 | 0.33 | -2188 | 2.12 | -2176 | 2.29 | -1987
MaxCut-D | $\\{0,1\\}$ | 0 | 0.5 | / | 0.24 | -22552 | 3.92 | -23904 | 0.34 | -20944 | 2.14 | -20616 | 2.28 | -16152
SK-Ising | $\\{-1,1\\}$ | 0 | 1 | complete | 0.21 | -44276 | 3.65 | -46280 | 0.34 | -43708 | 2.09 | -39904 | 2.24 | -31388
SK-QUBO | $\\{0,1\\}$ | 0 | 1 | complete | 0.18 | -22098 | 4.10 | -23982 | 0.33 | -22094 | 2.15 | -20266 | 2.09 | -15850
NAE-3SAT | $\mathbb{Z}$ | 0 | 0.0125 | / | 0.22 | -4012 | 0.95 | -4052 | 0.33 | -3388 | 2.11 | -3668 | 2.25 | -3268

We use the $1000$-spin $6$-regular Ising models for the performance
demonstration. The parameters of the Ising models are randomly generated
according to the uniform distribution
$\widehat{J}_{ij}\sim\mathcal{U}(-1,1),\widehat{h}_{i}\sim\mathcal{U}(-6,6)$,
Gaussian distribution
$\widehat{J}_{ij}\sim\mathcal{N}(0,\frac{1}{3}),\widehat{h}_{i}\sim\mathcal{N}(0,12)$,
and binary distribution
$\widehat{J}_{ij}\sim\mathcal{B}(\\{-\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}}\\}),\widehat{h}_{i}\sim\mathcal{B}(\\{-2\sqrt{3},2\sqrt{3}\\})$.
These three distributions have the same mean and variance. The experiment
results are shown in Tab. I, where we have repeated $1000$ trials for each
distribution. ADAM optimizer [39] is used with a learning rate of $6.0$ for
the gradient descent. The condition for convergence is set as

$\displaystyle\frac{\Delta\mathcal{L}(\bm{\theta})}{\mathcal{L}(\bm{\theta}_{\text{now}})-\mathcal{L}(\bm{\theta}_{\text{init}})}\leq
0.005,$ (51)

where $\bm{\theta}_{\text{now}},\bm{\theta}_{\text{init}}$ are the current and
initial $\bm{\theta}$. $\Delta\mathcal{L}(\bm{\theta})$ is the difference
between the maximum and the minimum of the losses in the last 5 iterations.
$\bm{\theta}_{\text{init}}$ is randomly initialized. For the uniform model,
the pruning rate is about $95\%$, which leads to a significant reduction in
computing time. The pruning also has a positive effect on the computing time
for Gaussian model. but the acceleration is smaller due to the smaller pruning
rate.

###  IV-C Comparisons

Here we use the time-to-solution (TTS) as a benchmark [40]

$\displaystyle T_{\text{TTS}}=T\frac{\log(1-0.99)}{\log(1-P)},$ (52)

where $P$ is the success probability of finding the optimal value in one run
and $T$ is the time cost of this run. Generally speaking, $P$ is difficult to
calculate since the optimal solution is unknown. Hence we use the estimation
$P=\frac{1}{d}e^{\chi}$ since the time complexity tends to be exponential
[41]. The estimated time-to-solution is then given by

$\displaystyle\widehat{T}_{\text{TTS}}=T\frac{\log(1-0.99)}{\log(1-\frac{1}{d}e^{\widehat{\chi}})}.$
(53)

Apart from TTS, the solution speed can also be simply defined as

$\displaystyle v=\frac{\widehat{\chi}}{T},$ (54)

which is taken as another metric of efficiency.

We compare our solver with the Python release of Gurobi v$10.0.1$, D-Wave
classical solver of steepest descent (DWave-SD), OpenJij with simulated
quantum annealing (OpenJij-SQA) and simulated annealing (OpenJij-SA). All the
experiments are conducted at the Apple M2 Pro core. To simulate the models
under broader situations, we generate the optimization problems using complete
graphs and the random parameters from the distributions
$J_{ij},h_{i}\sim\mathcal{U}(-1000,1000)$ and
$J_{ij},h_{i}\sim\mathcal{N}(0,\alpha^{2}),\alpha\sim\mathcal{U}(0,1000)$.
Comparison results are shown in Fig. 4, and more details can be found in Tab.
II.

(a) Comparisons on Key Indicators

(b) Statistical Visualization of Sample Points

Figure 4: Comparisons of different solvers on the randomly generated problems.
We sample $5000$ complete Ising problems of $1000$ spins. Except for Gurobi,
the number of sweeps is set to be $5$. For Gurobi, we set the gap of mixed-
integer programming as $\texttt{MIPGap}=0.8$. Running time is measured in
seconds.

From Fig. 4a and Tab. II we can see that our solver has the shortest solution
time in terms of TTS and $v$. In particular, TTS is a comprehensive metric
that characterizes the solution accuracy, the solution time, and the size of
the problem, which means that our solver our solver can find the global
optimal solution in the least time. The solution accuracy of our solver is the
second highest after Gurobi. Since the operation programming method
implemented in Gurobi often requires a large number of operations, its
computational complexity is high and the variance of its running time is
significantly higher than those of the other methods as can be seen in Tab.
II. Compared with the second fastest DWave-SD, our solver still has a clear
advantage. In the last, simulated annealing and quantum annealing do not show
any advantage over our solver in terms of performance, with worse solution
accuracy and a running time which is $9$ times more than that of ours.

To further illustrate the performance of our solver, we conduct experiments on
the typical benchmark optimization problems with the same settings, including
MaxCut-3 and MaxCut-D [42]. The Sherrington–Kirkpatrick (SK) model is an Ising
spin glass model with infinite spatial dimensions. We consider the SK spin-
glass model (SK-Ising) on a fully connected graph where the couplings
$J_{ij}=\pm 1$ are randomly chosen with equal probability, and SK-QUBO model
on a fully connected graph with the couplings $J_{ij}=0,1$ [41]. The
satisfiability problem (SAT) is one of the most fundamental NP-hard problems
which serves as a good benchmark for heuristic solvers. Not-All-Equal 3-SAT
(NAE-3SAT) is a variant of the Boolean SAT problem [42] for which we have used
randomly generated instances with the clause-to-variable ratio $2.11$ for
benchmarking. The experiment results can be found in Tab. III and Fig. 5.

Figure 5: Comparisons of different solvers on the typical fundamental
optimization problems.

First, we need to mention that pruning algorithms basically do not work on
these models with integer parameters. Our solver still holds the lead both in
accuracy and speed over the numerical and simulated quantum annealing methods.
Gurobi is more accurate than our solver. However, it can be seen that the
speed advantage of our solver over Gurobi increases as the problem densities
increase. For highly dense problems, the speed of our solver is generally $16$
times faster than that of Gurobi.

(a)

(b)

(c)

Figure 6: Comparisons of different solvers on scalability. Here we show the
variations of $T,\widehat{T}_{\text{TTS}},v$ with the number of spins $d$.
Apparently, the logarithmic horizontal and vertical axes can be related by
polynomial functions. The slope of the line is the highest power of the terms
of the polynomial, whose value is labeled on the right of the figure.

We have also explored the scalability of each solver by varying the size of
the Ising problems. The experiment results are summarized in Fig. 6, where we
have visualized the scalability in terms of the polynomial time complexity
with respect to $d$. It can be seen from Fig. 6a that the computational
complexity of our solver is only $\mathcal{O}(d^{0.59})$, which is much
smaller than the other schemes. In particular, Gurobi has the highest
computational complexity $\mathcal{O}(d^{2.18})$. The computational complexity
of the numerical methods and simulated quantum annealing is around
$\mathcal{O}(d^{2.0})$, suggesting that they have similar scalability
performance. DWave-SD is also significantly faster than the non-von Neumann
architectures.

##  V Conclusion

In this paper, we propose a numerical solver for the DIM. A pruning algorithm
has been designed which can effectively downsize the Ising model. An effective
domain selection method has been proposed to find the optimal solution in the
feasible domains by relaxing the discrete linear programming problem into an
optimization of parameters in a continuous space with gradient descent.

The experiments on the randomly generated and fundamental optimization
problems demonstrate that our solver is fast, efficient, and scalable. With
similar accuracy, our solver can be nearly $10$ times faster than the
classical Gurobi. Our scheme also shows a clear advantage in speed over the
steepest descent solver of D-Wave and the simulated annealing and simulated
quantum annealing algorithms of OpenJij. Therefore, our solver is more
suitable for speed-preferred applications, which has the potential to be
deployed in the near future to solve large-scale optimization problems, and in
the meantime serve as a preferred benchmark for evaluating the performance of
quantum computing.

## Acknowledgments

This work was supported by the National Key R&D Program of China (No.
2022YFB3304700) and the National Natural Science Foundation of China (No.
62173296).

## References

  * [1] E. Ising, “Contribution to the theory of ferromagnetism,” _Z. Phys_ , vol. 31, no. 1, pp. 253–258, 1925.
  * [2] J. Preskill, “Quantum computing in the nisq era and beyond,” _Quantum_ , vol. 2, p. 79, 2018. [Online]. Available: https://quantum-journal.org/papers/q-2018-08-06-79/?fbclid=IwAR0FRb9N2fas7ETWu2M40OS6prXB5QvFME_WRELpm2CAUcccVIEzA_UmLn4
  * [3] A. Lucas, “Ising formulations of many np problems,” _Frontiers in physics_ , vol. 2, p. 5, 2014. [Online]. Available: https://www.frontiersin.org/articles/10.3389/fphy.2014.00005/full?ref=https://githubhelp.com
  * [4] P. McMahon, “To solve optimization problems, just add lasers: An odd device known as an optical ising machine could untangle tricky logistics,” _IEEE Spectrum_ , vol. 55, no. 12, pp. 42–47, 2018. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/8544983
  * [5] S. J. Weinberg, F. Sanches, T. Ide, K. Kamiya, and R. Correll, “Supply chain logistics with quantum and classical annealing algorithms,” _Scientific Reports_ , vol. 13, no. 1, p. 4770, 2023. [Online]. Available: https://www.nature.com/articles/s41598-023-31765-8
  * [6] S. Bao, M. Tawada, S. Tanaka, and N. Togawa, “Multi-day travel planning using ising machines for real-world applications,” in _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_. IEEE, 2021, pp. 3704–3709. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/9564593
  * [7] Z. Pan, A. Sharma, J. Y.-C. Hu, Z. Liu, A. Li, H. Liu, M. Huang, and T. Geng, “Ising-traffic: Using ising machine learning to predict traffic congestion under uncertainty,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 37, no. 8, 2023, pp. 9354–9363. [Online]. Available: https://ojs.aaai.org/index.php/AAAI/article/view/26121
  * [8] H. Hussain, M. B. Javaid, F. S. Khan, A. Dalal, and A. Khalique, “Optimal control of traffic signals using quantum annealing,” _Quantum Information Processing_ , vol. 19, pp. 1–18, 2020. [Online]. Available: https://link.springer.com/article/10.1007/s11128-020-02815-1
  * [9] D. Inoue, A. Okada, T. Matsumori, K. Aihara, and H. Yoshida, “Traffic signal optimization on a square lattice with quantum annealing,” _Scientific reports_ , vol. 11, no. 1, p. 3303, 2021. [Online]. Available: https://www.nature.com/articles/s41598-021-82740-0
  * [10] D. Venturelli, S. Mandrà, S. Knysh, B. O’Gorman, R. Biswas, and V. Smelyanskiy, “Quantum optimization of fully connected spin glasses,” _Physical Review X_ , vol. 5, no. 3, p. 031040, 2015. [Online]. Available: https://journals.aps.org/prx/abstract/10.1103/PhysRevX.5.031040
  * [11] M. Parizy, P. Sadowski, and N. Togawa, “Cardinality constrained portfolio optimization on an ising machine,” in _2022 IEEE 35th International System-on-Chip Conference (SOCC)_. IEEE, 2022, pp. 1–6. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/9908082
  * [12] D. Venturelli and A. Kondratyev, “Reverse quantum annealing approach to portfolio optimization problems,” _Quantum Machine Intelligence_ , vol. 1, no. 1-2, pp. 17–30, 2019. [Online]. Available: https://link.springer.com/article/10.1007/s42484-019-00001-w
  * [13] M. Benedetti, J. Realpe-Gómez, R. Biswas, and A. Perdomo-Ortiz, “Estimation of effective temperatures in quantum annealers for sampling applications: A case study with possible applications in deep learning,” _Physical Review A_ , vol. 94, no. 2, p. 022308, 2016. [Online]. Available: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.94.022308
  * [14] F. W. Glover and G. A. Kochenberger, _Handbook of metaheuristics_. Springer Science & Business Media, 2006, vol. 57. [Online]. Available: https://link.springer.com/book/10.1007/b101874
  * [15] F. Böhm, D. Alonso-Urquijo, G. Verschaffelt, and G. Van der Sande, “Noise-injected analog ising machines enable ultrafast statistical sampling and machine learning,” _Nature Communications_ , vol. 13, no. 1, p. 5847, 2022. [Online]. Available: https://www.nature.com/articles/s41467-022-33441-3
  * [16] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, “A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem,” _Science_ , vol. 292, no. 5516, pp. 472–475, 2001. [Online]. Available: https://www.science.org/doi/abs/10.1126/science.1057726
  * [17] M. W. Johnson, M. H. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk _et al._ , “Quantum annealing with manufactured spins,” _Nature_ , vol. 473, no. 7346, pp. 194–198, 2011. [Online]. Available: https://www.nature.com/articles/nature10012
  * [18] T. Albash and D. A. Lidar, “Adiabatic quantum computation,” _Reviews of Modern Physics_ , vol. 90, no. 1, p. 015002, 2018. [Online]. Available: https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.90.015002
  * [19] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, “Coherent ising machine based on degenerate optical parametric oscillators,” _Physical Review A_ , vol. 88, no. 6, p. 063853, 2013. [Online]. Available: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.88.063853
  * [20] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu _et al._ , “A coherent ising machine for 2000-node optimization problems,” _Science_ , vol. 354, no. 6312, pp. 603–606, 2016. [Online]. Available: https://www.science.org/doi/abs/10.1126/science.aah4243
  * [21] V. Dunjko, J. M. Taylor, and H. J. Briegel, “Quantum-enhanced machine learning,” _Physical review letters_ , vol. 117, no. 13, p. 130501, 2016. [Online]. Available: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.117.130501
  * [22] M. Ohzeki, “Breaking limitation of quantum annealer in solving optimization problems under constraints,” _Scientific reports_ , vol. 10, no. 1, p. 3126, 2020. [Online]. Available: https://www.nature.com/articles/s41598-020-60022-5/
  * [23] S. Okada, M. Ohzeki, M. Terabe, and S. Taguchi, “Improving solutions by embedding larger subproblems in a d-wave quantum annealer,” _Scientific reports_ , vol. 9, no. 1, p. 2098, 2019. [Online]. Available: https://www.nature.com/articles/s41598-018-38388-4
  * [24] F. Böhm, T. Inagaki, K. Inaba, T. Honjo, K. Enbutsu, T. Umeki, R. Kasahara, and H. Takesue, “Understanding dynamics of coherent ising machines through simulation of large-scale 2d ising models,” _Nature communications_ , vol. 9, no. 1, p. 5020, 2018. [Online]. Available: https://www.nature.com/articles/s41467-018-07328-1
  * [25] S. Chowdhury, A. Grimaldi, N. A. Aadit, S. Niazi, M. Mohseni, S. Kanai, H. Ohno, S. Fukami, L. Theogarajan, G. Finocchio _et al._ , “A full-stack view of probabilistic computing with p-bits: devices, architectures and algorithms,” _IEEE Journal on Exploratory Solid-State Computational Devices and Circuits_ , 2023. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/10068500
  * [26] M. Yamaoka, C. Yoshimura, M. Hayashi, T. Okuyama, H. Aoki, and H. Mizuno, “A 20k-spin ising chip to solve combinatorial optimization problems with cmos annealing,” _IEEE Journal of Solid-State Circuits_ , vol. 51, no. 1, pp. 303–309, 2015. [Online]. Available: https://ieeexplore.ieee.org/document/7350099
  * [27] T. Takemoto, M. Hayashi, C. Yoshimura, and M. Yamaoka, “A 2 × 30k-spin multi-chip scalable cmos annealing processor based on a processing-in-memory approach for solving large-scale combinatorial optimization problems,” _IEEE Journal of Solid-State Circuits_ , no. 1, p. 55, 2020. [Online]. Available: https://ieeexplore.ieee.org/document/8906059
  * [28] K. Yamamoto, K. Kawamura, K. Ando, N. Mertig, T. Takemoto, M. Yamaoka, H. Teramoto, A. Sakai, S. Takamaeda-Yamazaki, and M. Motomura, “Statica: A 512-spin 0.25 m-weight annealing processor with an all-spin-updates-at-once architecture for combinatorial optimization with complete spin–spin interactions,” _IEEE Journal of Solid-State Circuits_ , vol. 56, no. 1, pp. 165–178, 2020. [Online]. Available: https://ieeexplore.ieee.org/document/9222223
  * [29] P. J. Van Laarhoven, E. H. Aarts, P. J. van Laarhoven, and E. H. Aarts, _Simulated annealing_. Springer, 1987. [Online]. Available: https://link.springer.com/book/10.1007/978-94-015-7744-1
  * [30] M. Yamaoka, C. Yoshimura, M. Hayashi, T. Okuyama, H. Aoki, and H. Mizuno, “24.3 20k-spin ising chip for combinational optimization problem with cmos annealing,” in _2015 IEEE International Solid-State Circuits Conference-(ISSCC) Digest of Technical Papers_. IEEE, 2015, pp. 1–3. [Online]. Available: https://ieeexplore.ieee.org/document/7063111
  * [31] N. A. Aadit, A. Grimaldi, M. Carpentieri, L. Theogarajan, J. M. Martinis, G. Finocchio, and K. Y. Camsari, “Massively parallel probabilistic computing with sparse ising machines,” _Nature Electronics_ , vol. 5, no. 7, pp. 460–468, 2022. [Online]. Available: https://www.nature.com/articles/s41928-022-00774-2
  * [32] Y. Haribara, S. Utsunomiya, and Y. Yamamoto, “A coherent ising machine for max-cut problems: performance evaluation against semidefinite programming and simulated annealing,” _Principles and Methods of Quantum Information Technologies_ , pp. 251–262, 2016. [Online]. Available: https://link.springer.com/chapter/10.1007/978-4-431-55756-2_12
  * [33] F. Rendl, G. Rinaldi, and A. Wiegele, “A branch and bound algorithm for max-cut based on combining semidefinite and polyhedral relaxations,” in _Integer Programming and Combinatorial Optimization: 12th International IPCO Conference, Ithaca, NY, USA, June 25-27, 2007. Proceedings 12_. Springer, 2007, pp. 295–309. [Online]. Available: https://link.springer.com/chapter/10.1007/978-3-540-72792-7_23
  * [34] M. T. Veszeli and G. Vattay, “Mean field approximation for solving qubo problems,” _Plos one_ , vol. 17, no. 8, p. e0273709, 2022. [Online]. Available: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0273709
  * [35] J. Bowles, A. Dauphin, P. Huembeli, J. Martinez, and A. Acín, “Quadratic unconstrained binary optimization via quantum-inspired annealing,” _Physical Review Applied_ , vol. 18, no. 3, p. 034016, 2022. [Online]. Available: https://journals.aps.org/prapplied/abstract/10.1103/PhysRevApplied.18.034016
  * [36] L. Li, H. Wang, Z. Xie, Z. Liu, W. Cui, and Y. Pan, “Simulated ising annealing algorithm with gaussian augmented hamiltonian monte carlo,” in _2023 42nd Chinese Control Conference (CCC)_. IEEE, 2023, pp. 6760–6765. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/10240679
  * [37] A. Bulut and T. K. Ralphs, “On the complexity of inverse mixed integer linear optimization,” _SIAM Journal on Optimization_ , vol. 31, no. 4, pp. 3014–3043, 2021. [Online]. Available: https://epubs.siam.org/doi/abs/10.1137/20M1377369
  * [38] N. L. Johnson, S. Kotz, and N. Balakrishnan, _Continuous univariate distributions, volume 2_. John wiley & sons, 1995, vol. 289. [Online]. Available: https://www.wiley.com/en-us/Continuous+Univariate+Distributions,+Volume+1,+2nd+Edition-p-9780471584957
  * [39] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014. [Online]. Available: https://arxiv.org/abs/1412.6980
  * [40] J. King, S. Yarkoni, M. M. Nevisi, J. P. Hilton, and C. C. McGeoch, “Benchmarking a quantum annealing processor with the time-to-target metric,” _arXiv preprint arXiv:1508.05087_ , 2015. [Online]. Available: https://arxiv.org/abs/1508.05087
  * [41] R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba, T. Honjo _et al._ , “Experimental investigation of performance differences between coherent ising machines and a quantum annealer,” _Science advances_ , vol. 5, no. 5, p. eaau0823, 2019. [Online]. Available: https://www.science.org/doi/10.1126/sciadv.aau0823
  * [42] H. Oshiyama and M. Ohzeki, “Benchmark of quantum-inspired heuristic solvers for quadratic unconstrained binary optimization,” _Scientific reports_ , vol. 12, no. 1, p. 2146, 2022. [Online]. Available: https://www.nature.com/articles/s41598-022-06070-5


# Stepped-Frequency THz-wave Signal Generation From a Kerr Microresonator
Soliton Comb

Omnia Nawwar, Kaoru Minoshima, and Naoya Kuse This work was supported by
Precursory Research for Embryonic Science and Technology (JPMJPR1905); Japan
Society for the Promotion of Science (21H01848, 21K18726); Cabinet Office,
Government of Japan (Subsidy for Reg. Univ. and Reg. Ind. Creation).O. Nawwar,
K. Minoshima, and N. Kuse are with Institute of Post-LED Photonics, Tokushima
University, Tokushima, Tokushima, 770-8506, Japan. K. Minoshima are also with
Graduate School of Informatics and Engineering, The University of Electro-
Communications, Chofu, Tokyo, 182-8585, Japan. (Corresponding author: Naoya
Kuse. email: kuse.naoya@tokushima-u.ac.jp).Manuscript received XXX; revised
XXX.

###### Abstract

Optically generated terahertz (THz) oscillators have garnered considerable
attention in recent years due to their potential for wide tunability and low
phase noise. Here, for the first time, a dissipative Kerr microresonator
soliton comb (DKS), which is inherently in a low noise state, is utilized to
produce a stepped-frequency THz signal ($\approx$ 280 GHz). The frequency of
one comb mode from a DKS is scanned through an optical-recirculating
frequency-shifting loop (ORFSL) which induces a predetermined frequency step
onto the carrier frequency. The scanned signal is subsequently heterodyned
with an adjacent comb mode, generating a THz signal in a frequency range that
is determined by the repetition frequency of the DKS. The proposed method is
proved by proof-of-concept experiments with MHz level electronics, showing a
bandwidth of 4.15 GHz with a frequency step of 83 MHz and a period of 16
$\mu$s.

###### Index Terms:

Linear stepped-frequency, Kerr microresonator frequency comb, terahertz
photonics, wideband waveform generation.

##  I Introduction

Figure 1: (a) Schematic of the experimental setup. AWG: arbitary waveform
generator, EDFA: Er-doped fiber amplifier, AOM: acousto-optic modulator, OBPF:
optical bandpass filter, UTC-PD: uni-travelling-carrier photodiode. (A) and
(B) show illustrations of the comb modes at locations (A) and (B). (C), (D)
and (E) show the instantaneous frequency of the comb mode/THz at locations
(C), (D) and (E). (b) Optical spectrum of the DKS. (c) Optical spectrum of the
filtered comb modes after the programmable OBPF.

The potential applications of terahertz (THz) signals (0.1-10 THz) have
experienced substantial expansion in recent years, fostering an escalating
interest in the advancement of THz technologies [1, 2, 3]. Specifically, the
necessity for wideband or frequency-scanned THz is unequivocal in applications
to radar [4]. Generation of THz signals can be principally accomplished via
two techniques - electronically or photonically. Electronic methods, wherein
microwave frequencies are multiplied to the THz domain, introduce inherent
obstacles due to nonlinearity, parasitic effects, and a decline in both system
efficiency and noise level, further complicated by fabrication difficulties
associated with achieving nanoscale precision. Alternatively, photonic
generation of THz signals has emerged as an advantageous substitute, where
optical frequency is downconverted into the THz domain. Photonic approaches
can be subdivided into two categories: photomixing of a single frequency
continuous-wave (CW) laser with a chirped-mode locked laser [5, 6], and
photomixing of two optical tones separated by the desired THz frequency [7, 8,
9, 10, 11, 12, 13, 14, 15]. As the former is encumbered by complexity and
limited reconfigurability (e.g., chirp rate), the latter has been extensively
adopted. For THz wave scanning, the frequency of one of the two optical tones
is varied. A simple strategy is to use two single-frequency CW lasers and
directly modulate the frequency of one [7, 8], but this transfers the CW
lasers’ phase noise to the generated THz wave, thereby degrading the
sensitivity of radar systems and hindering to reach scan-bandwidth-limited
resolution due to non-ideal linearity and reproducibility [16]. An alternative
approach to THz signal generation is through the use of two optical tones
extracted from electro-optic frequency combs [9, 10, 11, 12]. In this context,
the stability of the produced THz wave is only constrained by the reference
microwave signal, surpassing the performance of two independent CW lasers.
Nonetheless, the necessity for high-bandwidth electro-optic modulators (EOMs),
high-power microwave amplifiers, and a frequency-tunable microwave oscillator
to actuate the EOMs, convolute the entire system. In addition, the scan range
is restricted by the comb mode spacing (typically, 10 – 20 GHz) due to optical
bandpass filters (OBPFs) employed to isolate the two comb modes.

Recently, dissipative Kerr microresonator soliton combs (DKSs) [17, 18] have
been the focus of considerable attention as a method for THz wave generation
[13, 14, 15]. DKSs are generated by inputting a single-frequency CW laser into
high-Q microresonator, which are fabricated by CMOS-compatible processes [19],
thereby rendering them a chip-scale, mass-producible laser source [20].
Additionally, the comb mode spacing of DKSs ranges from 10 GHz – 1 THz, which
is suitable for the generation of THz waves. Furthermore, DKSs are in a mode-
locked state, showing high coherence among the comb modes. Along with locking
the repetition frequency of DKSs to external references such as a Brillouin
cavity [14] and fiber delay lines [21, 15], THz waves generated from DKSs
exhibits ultra-low phase noise with -100 dBc/Hz at a 10 kHz frequency offset
for 300 and 560 GHz carriers [14, 15]. However, the repetition frequency of
DKSs comb cannot be largely scanned. Even with a microheater deposited on a
micresonator, the scan range is 0.1 % of the repetition frequency [22, 23],
which corresponds to the frequency scanning of a 300 GHz wave of as small as
30 MHz, prohibiting the use of THz wave generated from DKSs for radar.

In this study, we propose and demonstrate a method to scan the frequency of a
THz wave generated from a DKS. Our proposed method extracts two neighboring
comb modes from a DKS and scans the frequency of one of these comb modes using
an optical recirculating frequency-shifting loop (ORFSL) [24, 25, 26, 27, 28,
29]. Upon heterodyning the two comb modes at a uni-traveling-carrier
photodiode (UTC-PD), a frequency-scanned THz wave is produced [30]. Our proof-
of-concept experiment effectively scans the frequency of a THz wave generated
from a DKS from 278.7 to 282. 8 GHz in 16 $\mu$s. This corresponds to a
bandwidth of 4.15 GHz with a frequency step of 83 MHz.

##  II Experimental setup and operation principle

A basic architecture of the experimental setup and signal at different points
are depicted in Fig. 1(a). The output of a pump CW laser is modulated by a
dual-parallel Mach-Zehnder modulator (DP-MZM) (not shown in Fig. 1(a)),
amplified by an Er-doped fiber amplifier (EDFA) (not shown in Fig. 1(a)), and
coupled into a high-Q Si3N4 microresonator (Ligentec SA) with a free-spectral
range of about 280 GHz. The DP-MZM is operated in a carrier-suppressed single-
sideband (CS-SSB) mode, which is used to rapidly scan the frequency of the
pump CW laser to access a stable DKS [31]. More details to generate a DKS is
described in the appendix and ref [32]. The optical spectrum of the DKS used
in this work is shown in Fig. 1(b). The comb mode spacing and 10-dB bandwidth
are about 280 GHz and 90 nm, respectively.

A bandstop filter (not shown in Fig. 1(a)) is used to reject the residual pump
CW laser before feeding the DKS to a programmable OBPF to pass a pair of
neighboring comb modes used in our experiment. The two neighboring comb modes
at the wavelengths of 1557.44 and 1555.14 nm, corresponding to -6th and -5th
with respect to the pump mode, are shown in Fig. 1(c). The signal-to-noise
ratio (SNR) of the comb modes is about 38 dB with the resolution bandwidth
(RBW) of 0.02 nm. Extracted modes are amplified to about 60 mW before
splitting into two branches through a 50/50 optical coupler. Individual comb
modes are extracted by OBPFs in the upper and lower branches as shown in Fig.
1(a) at locations (A) and (B). The signal in the upper branch experiences the
frequency shifts by an ORFSL [33]. The ORFSL consists of an acousto-optic
modulator (AOM) (AOM 1 in Fig. 1(a)) with an extinction ratio of more than 50
dB and a frequency shifting loop (FSL). AOM 1, which is drived by a RF signal
from an arbitrary waveform generator (AWG), works as an optical switch to
convert the comb mode into optical pulses. The pulse width equals to the time
delay provided by the ORFSL, and the pulse cycle is determined by the number
of round trips (N) of the ORFSL.

The generated pulses are directed into the ORFSL which is realized by a 50/50
optical coupler, fiber delay, EDFA 2, OBPF, AOM 2, and polarization
controller. Apart from a 8 cm fiber on the polarization controller, fibers and
fiber components in the ORFSL are polarization maintaining. The optical
coupler is used as input and output ports of the ORFSL. The length of the
fiber used is 49 m which determines the time delay in the loop. EDFA 2
compensates the loss in the loop (mainly due to coupling and insertion losses
for components) while the following BPF suppresses the amplified spontaneous
emission (ASE) noise. AOM 2 shifts the pulse frequency for every round trip by
a frequency step of $\Delta f$ (= 83 MHz in our experiment), which corresponds
to the driving frequency from the AWG, ensuring the time-frequency linearity
of resultant stepped-frequency THz signal. When AOM 2 is activated, the input
pulse successively experiences the frequency shift of $\Delta f$ every round
trips. AOM 2 is deactivated when a new incoming pulse from AOM 1 enters the
ORFSL, thereby initiating a fresh cycle of stepped-frequency signal generation
process. Note that the driving signals of AOMs 1 and 2 are synchronized so
that the timing of the operation of AOMs 1 and 2 are precisely controlled. For
each cycle, the optical pulse experiences a frequency shift of $N\times\Delta
f$ as shown in Fig. 1(a) at location (C). Another 50/50 optical coupler is
used to combine the scanned signal from the upper branch and the comb mode
from the lower branch, whose instantaneous frequency is constant at location
(D). Then, it is further amplified by an EDFA (not shown in Fig. 1(a) to
provide about 20 mw optical power to UTC-PD (IOD-PMJ-13001, 280 - 380 GHz).

Owing to the square-law envelope detection at the UTC-PD, the UTC-PD generates
a frequency-stepped THz signal by down-converting the optical frequencies of
the comb modes to a THz wave as shown in Fig. 1(a) at location (E). The
frequency of the THz wave corresponds to the spacing between the scanned comb
mode and the neighboring comb mode, resulting in the scan range of
$N\times\Delta f$. For the following experimental results, the time delay and
frequency step are kept the same.

##  III Experimental results

First, we examine a single tone THz signal generated from either two comb
modes or two free-running CW lasers without implementing the ORFSL. Our
investigation primarily focuses on the SNR and power of the THz signals. The
produced THz signal is characterized by down-converting the THz signal to a
microwave using the setup shown in Fig. 2(a). The THz from the UTC-PD is mixed
with a frequency-multiplied ($\times$ 12) local oscillator (11.79 GHz)
followed by an amplifier (the frequency multiplier and amplifier are
incorporated into a single device, WR6.5AMC-I from Virginia Diodes, Inc.) at a
sub-harmonic mixer (SHM, WR3.4SHM from Virginia Diodes, Inc.). The down-
converted signal is then further amplified by microwave amplifiers and
digitally detected by an electrical spectrum analyzer (ESA). Figure 2 (b)
shows the relationship between the optical power to the UTC-PD and RF power
after down conversion. The THz power (proportional to RF power) is shown to be
proportional to the square of the optical input power, mirroring the principle
of standard square-law envelope PDs. Although we don’t measure the THz power
directly, it is estimated to be 10 $\mu$W according to the datasheet, which
could potentially be increased up to 100 $\mu$W by augmenting the optical
input power. There is a little difference in power between the down-converted
RF signals from the two comb modes and two free-running CW lasers. There is a
minor difference in power between the down-converted RF signals originating
from the two comb modes and the two free-running CW lasers. This discrepancy
is not fundamental, but likely stems from measurement errors or system
reproducibility issues such as polarization and temperature dependence of the
UTC-PD response. Indeed, the photocurrent of the UTC-PD shows a slight
variation from day to day, even with consistent optical input power. The RF
spectra of the signals from two CW lasers and comb modes are depicted in Figs.
2(c) and (d). Both cases exhibit almost the same power and noise floor,
suggesting that despite the optical SNR of the comb modes being lower than
that of the two CW lasers, the quality (i.e., power and SNR) of the
photonically generated THz signal is not compromised by using the comb modes.
Currently, the noise floor of -65 dBm is limited by the leakage of the local
oscillator from the SHM.

Figure 2: (a) Schematic of the system to down convert a THz signal to a
microwave. SHM: sub-harmonic mixer, ESA: electrical spectrum analyzer. (b)
Power of down-converted signal when two CW lasers (blue triangle) and two comb
modes (red circle) are used. The dotted line show a fitting with a square
function to the result of the power when two comb modes are used. (c) RF
spectrum of the down-converted signal (blue curve) and noise floor (black
curve) when two CW lasers are used. (d) RF spectrum of the down-converted
signal (blue curve) and noise floor (black curve) when two comb modes are
used.

Figure 3: (a) (Blue curve) The intensity of the outputs from the ORFSL. (Red
curve) The down-converted signal from a generated THz when the number of round
trips in the ORFSL is 32. One cycle is highlighted by the blue square.
(b)(c)(d) FFT of the down-converted signal at 1st, 16th, and 32nd,
respectively, showing the instanteneous frequency of 0.165 GHz, 1.41 GHz, and
2.738 GHz, respectively. (e) Spectrogram of the down-converted signal when the
number of round trips in the ORFSL is 50.

Next, we turn our attention to the generation of step-frequency signals
utilizing two comb modes. The comb modes employed in this process are
extracted by the programmable OBPF as depicted in Fig. 1(c). The mode with a
higher frequency is channeled into the ORFSL with the period of the RF signal
to the AOMs set to 10.24 µs (equivalent to 32 round trips) and a pulse width
of 320 ns. This pulse width is intended to match the time delay within the
ORFSL and is experimentally adjusted to minimize the intensity disparity
across time slots. The intensity of the output from the ORFSL is measured by
subtly coupling out the light after the ORFSL (not shown in Fig. 1(a)), which
is directed to a PD (not shown in Fig. 1(a)). Blue curve in Fig. 3(a)
represents the optical pulses detected after the ORFSL with the PD while the
red curve depicts the heterodyned signal detected with the UTC-PD and
subsequently down-converted to be displayed on a fast oscilloscope. Given the
frequency’s gradual increase from the initial to the final loop, later pulses
undergo greater attenuation than their leading counterparts due to the
frequency-dependent responses of THz components. To counterbalance this
attenuation, the intensity of the optical pulses preceding the UTC-PD is
escalated, as measured by the PD and displayed as the blue curve in Fig. 3(a).
This adjustment is achieved by modulating the pump current of EDFA 2 and the
polarization controller within the ORFSL. Figures 3(b), (c), and (d)
illustrate the 1st, 16th, and 32nd frequency-shifting pulses within the time
domain, displaying a clear frequency increase. The frequency of these bursts
(post-down conversion) is calculated using a Fast Fourier Transform (FFT) to
be 0.165 GHz, 1.41 GHz, and 2.738 GHz, respectively, each with a frequency
step of 83 MHz, thereby yielding a signal with a bandwidth of 2.537 GHz.
Figure 4 presents a spectrogram showcasing the instantaneous frequency of all
pulses when a total of 50 round trips (corresponding to a period of 16 $\mu$s)
are utilized. The dwell time within one time slot remains at 320 ns, as
dictated by the time delay within the loop. The instantaneous frequency of the
steps consistently increases with a frequency step of 83 MHz from 0.166 MHz to
4.233 GHz. When considering the frequency down-conversion, the THz frequency
is scanned from 278.727 GHz to 282.793 GHz, achieving a bandwidth of 4.1 GHz.

Lastly, we explore the feasibility of accommodating additional round trips.
Two primary concerns arise. First, there’s the potential escalation of the
system’s noise floor, which could be induced by the accumulated noise from the
EDFA, coupled with the minimal suppression of the unwanted carrier from the
AOM inside the loop. Second, there’s the risk of amplitude fluctuations,
likely due to the accumulation of imperfect polarization maintaining
procedures, paired with the presence of polarization-dependent components in
the loop, such as an isolator and Er fiber. Figure 4 shows FFT signals of a
single time slot at the 1st, 16th, and 32nd rounds. It is noteworthy that the
noise floor does not escalate with an increased number of round trips. In this
context, the RF power is not a concern as it is solely determined by the
optical input power to the UTC-PD. This outcome could be specific to the THz
system, where the noise floor is constrained by the detection setup rather
than by the optical tones. Indeed, the degradation of the SNR has been
observed in the system for the stepped-frequency microwave using an ORFSL, as
mentioned in reference [24, 26].

Figure 4: RF spectra of FFT signals of a down-converted signal at 1st, 16th,
and 32nd time slots when the number of round trips is 32.

While it would be beneficial to observe the noise floor across a larger number
of round trips, the limited bandwidth of the equipment (such as the
oscilloscope and microwave amplifier) available in our lab precludes
conducting the same experiments with more round trips. However, we can assess
amplitude fluctuation without generating THz waves by measuring the intensity
of the output from the ORFSL. Optical pulses for 100 and 400 round trips are
illustrated in Figs. 5(a) and (b), respectively. With a total of 100 round
trips, the intensity can be equalized by carefully adjusting the pump current
of the EDFA and the polarization controller within the ORFSL. However,
managing the optical output power from the loop becomes exceedingly
challenging when the total number of round trips reaches 400. The scan
bandwidth of the THz is anticipated to be 8.3 GHz for 100 round trips, and
33.2 GHz for 400 round trips, respectively.

Figure 5: (a) and (b) Intensity of the outputs from the ORFSL when the number
of round trips are 100 and 400, respectively. One cycle is highlighted by the
blue square.

##  IV Discussions and conclusion

While we maintain a fixed time delay and frequency step in the experiments,
both parameters can be adjusted. The frequency step of the ORFSL, for
instance, can be altered from several MHz to tens of GHz. This can be achieved
by choosing an AOM with a different modulation frequency, cascading multiple
AOMs, or using a DP-MZM operating in the carrier-suppressed single-sideband
mode [27]. As a result, the ORFSL offers a versatile frequency step
capability, adaptable to meet the specific spectral resolution requirements of
various applications. Similarly, the dwell time of the step-frequency THz
signal within a single time slot can be adjusted by altering the length of the
fiber loop, which consequently controls the pulse width. A pulse width of a
few nanoseconds can be introduced to achieve shorter dwell times. Notably,
since the frequency step and dwell time are governed by separate experimental
components, these parameters can be independently optimized, offering a high
degree of design flexibility. Furthermore, the number of round trips can be
controlled by simply adjusting the driving RF signals to the AOMs, resulting
in changes to the bandwidth of the step-frequency THz signals. However, such
escalations are not unlimited. Theoretically, while power loss incurred during
round trips can be offset by intra-loop EDFA, the SNR continuously decreases
until it reaches an unacceptable level. Nevertheless, the SNR of the two modes
from the DKSs does not limit the quality of the generated THz, at least up to
the 32nd round trip. The SNR of the comb modes can be improved by utilizing
the injection locking of DFB lasers to the comb modes, without sacrificing
phase noise and frequency stability [34]. More round trips also lead to
amplitude fluctuations, as observed in Fig. 5(b). One strategy to counteract
these fluctuations involves installing feedback loops to stabilize pulse
energy from the loop by controlling the modulation amplitude of AOM 2 and the
pump current of EDFA 2. Alternatively, the observed amplitude fluctuations can
be used to normalize the THz power by simply dividing by the square of the
amplitude fluctuation, though in practice, additional minor modification
factors may be necessary.

In conclusion, we successfully generated a step-frequency-THz signal from
DKSs, in which one comb mode’s frequency was discretely scanned using an
ORFSL, and then heterodyned with a neighboring comb mode. Our experiment
involved an ORFSL with 32 and 50 frequency steps (one frequency step = 83
MHz), generating a THz signal with a start frequency of 278.7 GHz and a
bandwidth of 4.1 GHz. We also investigated the feasibility of increasing the
number of ORFSL round trips, with a particular focus on the SNR and amplitude
fluctuation. DKSs demonstrate a high degree of phase correlation between the
comb modes, resulting in reduced phase noise distortion and frequency
fluctuation of THz signals. The developed system holds potential for
application in THz radar [35, 24], which could be pivotal technology in the
era of beyond 5G or 6G.

##  Appendix A Detailed experimental setup

###  A-A DKS generation

A pump CW laser operating at wavelength of 1543 nm is utilized. To rapidly
scan the frequency of the pump CW laser, the CW laser passes through a DP-MZM.
The DP-MZM is modulated by a voltage-controlled oscillator (VCO). The output
of the VCO is amplified and then divided into two signals with a 90-degree
phase difference using a 90-degree splitter. These split RF signals are
subsequently applied to the DP-MZM. Once the DC biases for the DP-MZM are
properly adjusted, the DP-MZM operates in the CS-SSB mode. By altering the
frequency of the VCO, the DP-MZM functions as a frequency shifter. The
frequency of the VCO is changed by a few GHz from 10 GHz in less than 100 ns.
The output from the DP-MZM is then amplified using an EDFA, followed by an
OBPF to eliminate ASE from the EDFA. Finally, the pump CW laser, with a power
level of about 300 mW, is coupled into a chip with the a microresonator,
experiencing a coupling loss of `<` 3 dB.

###  A-B Clock signal generation

The RF signals applied to the AOMs are produced by AWGs and RF switches. The
AWG generates two outputs: a sine wave with a frequency suitable for the AOMs
(approximately 80 MHz in our case) and a pulse train. The pulse width and
period are determined by the time delay and the number of round trips in the
FSL. These sine wave and pulse train are fed into an RF switch, which
generates a pulse train with a carrier frequency that corresponds to the
frequeny of the sine wave.

Figure 6: Schematic of the experimental setup. VCO: voltage-controlled
oscillator.

##  Appendix B Micro/mm/THz wave generation using ORFSL

The main claim of this study is the generation of frequency-stepped THz signal
from a DKS. However, we believe it is instructive to revisit the literature on
the generation of micro/mm/THz waves generation using ORFSLs. Table I presents
a summary of the essential parameters. For microwave (`<` 30 GHz) generation
[26, 28, 29, 27], a single cw laser is utilized. The CW laser is split, and
the frequency of one of the two resultant beams is shifted by the ORFSL, which
enables scanning from approximately DC upto the bandwidth. The phase noise of
the CW laser is cancelled out when a microwave is generated. For mm wave
($\approx$ 30 GHz) [24], two CW laser with a frequency offset are employed to
add a frequency offset to the generated mm wave. Due to the use of two
independent CW lasers, the phase noise of these lasers is transferred to the
generated mm wave. For THz generation ($\approx$ 300 GHz) as demonstrated in
this work, a significant frequency offset is added by using two optical tones
from a DKS. Since the DKS operates in a mode-locked state, the phase noise of
the THz wave (equivalent to relative phase noise between the comb modes) is
far superior to the case where a THz generated from two CW lasers.

TABLE I: KEY PARAMETERS OF THE LINEAR FREQUENCY STEPPED SIGNALS Ref | Source of two tones | Center Frequency (GHz) | BW (GHz) | Step frequency$\Delta f$ (GHz) | Max. N
---|---|---|---|---|---
[26] | One CW laser | 10 | 20 | 0.04 | 500
[28] | One CW laser | 12.5 | 25 | 0.085 | 300
[29] | One CW laser | 15 | 30 | 0.1 | 200
[27] | One CW laser | 20 | 26 | 1, 3, 8, 10 | 29
[24] | Two CW laser | 34 | 20 | 0.04/0.08 | 144
This work | DKS | 281 | 4.1 | 0.083 | 50

## References

  * [1] P. A. Banks, E. M. Kleist, and M. T. Ruggiero, “Investigating the function and design of molecular materials through terahertz vibrational spectroscopy,” Nat. ReviewsChemistry (2023).
  * [2] X. Chen, H. Lindley-Hatcher, R. I. Stantchev, J. Wang, K. Li, A. H. Serrano, Z. D. Taylor, E. Castro-Camus, and E. Pickwell-MacPherson, “Terahertz (THz) biophotonics technology: Instrumentation, techniques, and biomedical applications,” Chemical Phys.Reviews 3, 011311 (2022).
  * [3] T. Nagatsuma, G. Ducournau, and C. C. Renaud, “Advances in terahertz communications accelerated by photonics,” Nat.Photonics 10, 371–379 (2016).
  * [4] K. B. Cooper, R. J. Dengler, N. Llombart, B. Thomas, G. Chattopadhyay, and P. H. Siegel, “THz imaging radar for standoff personnel screening,” IEEE Transactions on Terahertz Science Technology 1, 169–182 (2011).
  * [5] I. Lin, J. McKinney, and A. Weiner, “Photonic synthesis of broadband microwave arbitrary waveforms applicable to ultra-wideband communication,” IEEE Microw. Wirel. ComponentsLetters 15, 226–228 (2005).
  * [6] A. Rashidinejad and A. M. Weiner, “Photonic radio-frequency arbitrary waveform generation with maximal time-bandwidth product capability,” J. LightwaveTechnology 32, 3383–3393 (2014).
  * [7] G. Danion, C. Hamel, L. Frein, F. Bondu, G. Loas, and M. Alouini, “Dual frequency laser with two continuously and widely tunable frequencies for optical referencing of GHz to THz beatnotes,” OpticsExpress 22, 17673 (2014).
  * [8] E. A. Kittlaus, D. Eliyahu, S. Ganji, S. Williams, A. B. Matsko, K. B. Cooper, and S. Forouhar, “A low-noise photonic heterodyne synthesizer and its application to millimeter-wave radar,” Nat.Communications 12 (2021).
  * [9] A. Hirata, H. Togo, N. Shimizu, H. Takahashi, K. Okamoto, and T. Nagatsuma, “Low-phase noise photonic millimeter-wave generator using an awg integrated with a 3-dB combiner,” IEICE transactions onelectronics 88, 1458–1464 (2005).
  * [10] B. Jerez, F. Walla, A. Betancur, P. Martín-Mateos, C. de Dios, and P. Acedo, “Electro-optic THz dual-comb architecture for high-resolution, absolute spectroscopy,” OpticsLetters 44, 415 (2019).
  * [11] A. Parriaux, K. Hammani, and G. Millot, “Electro-optic frequency combs,” Advances Optics Photonics 12, 223–287 (2020).
  * [12] S. Jia, M.-C. Lo, L. Zhang, O. Ozolins, A. Udalcovs, D. Kong, X. Pang, R. Guzman, X. Yu, S. Xiao, S. Popov, J. Chen, G. Carpintero, T. Morioka, H. Hu, and L. K. Oxenløwe, “Integrated dual-laser photonic chip for high-purity carrier generation enabling ultrafast terahertz wireless communications,” Nat.Communications 13 (2022).
  * [13] S. Zhang, J. M. Silver, X. Shang, L. D. Bino, N. M. Ridler, and P. Del’Haye, “Terahertz wave generation using a soliton microcomb,” OpticsExpress 27, 35257 (2019).
  * [14] T. Tetsumoto, T. Nagatsuma, M. E. Fermann, G. Navickaite, M. Geiselmann, and A. Rolland, “Optically referenced 300 GHz millimetre-wave oscillator,” Nat.Photonics 15, 516–522 (2021).
  * [15] N. Kuse, K. Nishimoto, Y. Tokizane, S. Okada, G. Navickaite, M. Geiselmann, K. Minoshima, and T. Yasui, “Low phase noise THz generation from a fiber-referenced Kerr microresonator soliton comb,” CommunicationsPhysics 5 (2022).
  * [16] T.-J. Ahn, J. Y. Lee, and D. Y. Kim, “Suppression of nonlinear frequency sweep in an optical frequency-domain reflectometer by use of Hilbert transformation,” Appl.optics 44, 7630–7634 (2005).
  * [17] T. Herr, V. Brasch, J. D. Jost, C. Y. Wang, N. M. Kondratiev, M. L. Gorodetsky, and T. J. Kippenberg, “Temporal solitons in optical microresonators,” Nat.Photonics 8, 145–152 (2014).
  * [18] T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, “Dissipative Kerr solitons in optical microresonators,” Science 361, eaan8083 (2018).
  * [19] A. Kovach, D. Chen, J. He, H. Choi, A. H. Dogan, M. Ghasemkhani, H. Taheri, and A. M. Armani, “Emerging material systems for integrated optical Kerr frequency combs,” Advances Optics Photonics 12, 135–222 (2020).
  * [20] C. Xiang, J. Liu, J. Guo, L. Chang, R. N. Wang, W. Weng, J. Peters, W. Xie, Z. Zhang, J. Riemensberger, J. Selvidge, T. J. Kippenberg, and J. E. Bowers, “Laser soliton microcombs heterogeneously integrated on silicon,” Science 373, 99–103 (2021).
  * [21] D. Kwon, D. Jeong, I. Jeon, H. Lee, and J. Kim, “Ultrastable microwave and soliton-pulse generation from fibre-photonic-stabilized microcombs,” Nat.Communications 13, 1–8 (2022).
  * [22] N. Kuse, T. Tetsumoto, G. Navickaite, M. Geiselmann, and M. E. Fermann, “Continuous scanning of a dissipative Kerr-microresonator soliton comb for broadband, high-resolution spectroscopy,” OpticsLetters 45, 927 (2020).
  * [23] N. Kuse, G. Navickaite, M. Geiselmann, T. Yasui, and K. Minoshima, “Frequency-scanned microresonator soliton comb with tracking of the frequency of all comb modes,” OpticsLetters 46, 3400–3403 (2021).
  * [24] Y. Liu, Z. Zhang, M. Burla, and B. J. Eggleton, “11-GHz-bandwidth photonic radar using MHz electronics,” Laser $\&$ PhotonicsReviews 16, 2100549 (2022).
  * [25] T. Chen, W. Kong, H. Liu, and R. Shu, “Frequency-stepped pulse train generation in an amplified frequency-shifted loop for oxygen a-band spectroscopy,” OpticsExpress 26, 34753 (2018).
  * [26] Y. Lyu, Y. Li, C. Yu, L. Yi, T. Nagatsuma, and Z. Zheng, “Photonic generation of highly-linear ultra-wideband stepped-frequency microwave signals with up to 6.106 time-bandwidth product,” J. LightwaveTechnology 40, 1036–1042 (2022).
  * [27] Y. Zhang, C. Liu, K. Shao, Z. Li, and S. Pan, “Multioctave and reconfigurable frequency-stepped radar waveform generation based on an optical frequency shifting loop,” OpticsLetters 45, 2038 (2020).
  * [28] C. Schnébelin, J. Azaña, and H. G. de Chatellus, “Programmable broadband optical field spectral shaping with megahertz resolution using a simple frequency shifting loop,” Nat.Communications 10 (2019).
  * [29] Z. Zhang, Y. Liu, and B. J. Eggleton, “Photonic generation of 30 GHz bandwidth stepped-frequency signals for radar applications,” J. LightwaveTechnology 40, 4521–4527 (2022).
  * [30] H. Ito, T. Furuta, F. Nakajima, K. Yoshino, and T. Ishibashi, “Photonic generation of continuous THz wave using uni-traveling-carrier photodiode,” J. LightwaveTechnology 23, 4016–4021 (2005).
  * [31] J. R. Stone, T. C. Briles, T. E. Drake, D. T. Spencer, D. R. Carlson, S. A. Diddams, and S. B. Papp, “Thermal and nonlinear dissipative-soliton dynamics in Kerr-microresonator frequency combs,” Physical reviewletters 121, 063902 (2018).
  * [32] N. Kuse, T. C. Briles, S. B. Papp, and M. E. Fermann, “Control of Kerr-microresonator optical frequency comb by a dual-parallel Mach-Zehnder interferometer,” OpticsExpress 27, 3873–3883 (2019).
  * [33] K. Shimizu, T. Horiguchi, and Y. Koyamada, “Technique for translating light-wave frequency by using an optical ring circuit containing a frequency shifter,” Opticsletters 17, 1307–1309 (1992).
  * [34] N. Kuse and K. Minoshima, “Amplification and phase noise transfer of a Kerr microresonator soliton comb for low phase noise THz generation with a high signal-to-noise ratio,” OpticsExpress 30, 318 (2021).
  * [35] C. Wang, Q. Zhang, J. Hu, C. Li, S. Shi, and G. Fang, “An efficient algorithm based on CSA for THz stepped-frequency SAR imaging,” IEEE Geoscience Remote Sens.Letters 19, 1–5 (2022).