Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.
10.1109/ACCESS.2023.3310875

Corresponding author: Shuta Kikuchi (e-mail<EMAIL_ADDRESS>
# Dynamical process of a bit-width reduced Ising model with simulated
annealing

SHUTA KIKUCHI1  NOZOMU TOGAWA2  and SHU TANAKA1,3,4,5  Department of Applied
Physics and Physico-Informatics, Keio University, Kanagawa 223-8522, Japan
Department of Computer Science and Communications Engineering, Waseda
University, Tokyo 169-8555, Japan  Human Biology-Microbiome-Quantum Research
Center (WPI-Bio2Q), Keio University, Tokyo 108-8345, Japan  Green Computing
System Research Organization, Waseda University, Tokyo 162-0042, Japan
International Research Frontiers Initiative, Tokyo Institute of Technology,
Tokyo, 108-0023, Japan

###### Abstract

Ising machines have attracted attention as efficient solvers for combinatorial
optimization problems, which are formulated as ground-state (lowest-energy)
search problems of the Ising model. Due to the limited bit-width of
coefficients on Ising machines, the Ising model must be transformed into a
bit-width reduced (BWR) Ising model. According to previous research, the bit-
width reduction method, which adds auxiliary spins, ensures that the ground
state of the BWR Ising model is theoretically the same as the Ising model
before bit-width reduction (original Ising model). However, while the
dynamical process is closely related to solution accuracy, how the BWR Ising
model progresses towards the ground state remains to be elucidated. Therefore,
we compared the dynamical processes of these models using simulated annealing
(SA). Our findings reveal significant differences in the dynamical process
across models. Analysis from the viewpoint of statistical mechanics found that
the BWR Ising model has two characteristic properties: an effective
temperature and a slow relaxation. These properties alter the temperature
schedule and spin flip probability in the BWR Ising model, leading to
differences in the dynamical process. Therefore, to obtain the same dynamical
process as the original Ising model, we proposed SA parameters for the BWR
Ising model. We demonstrated the proposed SA parameters using a square lattice
Ising model, in which all coefficients were set uniformly to the same positive
values or randomly. Our experimental evaluations demonstrated that the
dynamical process of the BWR and original Ising model became closer.

###### Index Terms:

Bit-width reduction, Ising machine, Ising model, simulated annealing,
statistical mechanics

=-21pt

##  I Introduction

###  I-A Combinatorial optimization problem and Ising model

Combinatorial optimization problems find the optimal combination of decision
variables to minimize or maximize the objective function for the given
constraints. Typical examples include the traveling salesman problem, the Max-
Cut problem, and the knapsack problem. Because such problems can be found in
many real-world application domains, there is growing interest in developing
techniques to find the optimal or quasi-optimal solution efficiently and
accurately.

Some combinatorial optimization problems can be formulated in a mathematically
constructed model in statistical mechanics called an Ising model or its
equivalent model called a quadratic unconstrained binary optimization (QUBO)
model [1, 2]. The ground state of the Ising model corresponds to the optimal
solution of the combinatorial optimization problem, where the ground state is
referred to as the lowest-energy state.

An Ising model is defined on an undirected graph $G=(V,E)$, where $V$ and $E$
are sets of vertices and edges, respectively. The Ising model consists of
spins, magnetic fields, and interactions. The Hamiltonian (or energy function)
$H$ of the Ising model is defined by

$\displaystyle H=-\sum_{i\in V}h_{i}\sigma_{i}-\sum_{(i,j)\in
E}J_{ij}\sigma_{i}\sigma_{j},$ (1)

where $\sigma_{i}$ is the spin on the vertex $i\in V$ and has a value of
$(+1)$ or $(-1)$. $h_{i}$ is the magnetic field on the vertex $i\in V$, and
$J_{ij}$ is the interaction on the edge $(i,j)\in E$. In this paper, we assume
that interactions and magnetic fields are integer constants.

###  I-B Ising machine

Approaches such as meta-heuristics and Ising machines have been developed to
solve combinatorial optimization problems [3, 4, 5, 6, 7, 8, 9, 10, 11, 12].
Ising machines have attracted attention as fast and high-precision solvers for
combinatorial optimization problems. Ising machines specialize in searching
for better solutions to combinatorial optimization problems formulated by an
Ising model or a QUBO model. Studies have applied Ising machines to various
combinatorial optimization problems, including machine learning [13, 14, 15,
16], material design [17, 18, 19], portfolio optimization [20, 21], protein
folding [22], traffic optimization [23, 24, 25, 26, 27], quantum compiler
[28], and black-box optimization [17, 29, 30].

For an Ising machine to solve the problems formulated in the Ising model, the
model must be mapped to the machine [31]. However, Ising machines are limited
by their hardware specifications. For example, the number of spins
corresponding to the problem size, the topology related to the connectivity
between spins, and the bit-width which is the imputable numerical range for
coefficients of the interactions $J_{ij}$ and magnetic fields $h_{i}$. The
specifications of the various Ising machines are summarized in [32, 33].
Although various approaches have been proposed in previous research to
overcome the limitations due to the number of spins [34, 35, 36, 37, 38, 39,
40] and the topology [41, 42, 43, 44, 45, 46, 47], few studies have been
devoted to overcoming the bit-width limitation.

###  I-C Motivation of this study

The bit-width of a digital Ising machine, implemented by digital circuits such
as Graphics Processing Unit (GPU), Field Programmable Gate Array (FPGA), or
Application Specific Integrated Circuit (ASIC), is represented by an integer
range of sign bits. Here, we assume that a bit-width of $n$-bits shows
$[-(2^{n-1}-1),2^{n-1}-1]$.

When the coefficients of the Ising model exceed the implemented bit-width of
the Ising machine, they cannot be inputted into the Ising machine. Thus, bit-
width reduction methods such as the shift method are used. The shift method
divides by two until the coefficients of the Ising model fall within the
target bit-width range [32]. Although the shift method can naively reduce bit-
width, the ground states of the bit-width reduced (BWR) Ising model may differ
from that of the Ising model before bit-width reduction (original Ising
model). Therefore, a new method to reduce bit-width by adding auxiliary spins
is proposed [32]. The proposed method guarantees that the ground states of the
original Ising model and the BWR Ising model are theoretically consistent.
However, the dynamical process of the BWR Ising model towards the ground state
remains to be elucidated.

This study analyzes the dynamical process of the BWR Ising model using
simulated annealing (SA), which is the most fundamental algorithm for Ising
machines implemented with digital circuits. The contributions of this study
are as follows:

  * •

The difference between the dynamical processes of the original Ising model and
that of the BWR Ising model applying the proposed bit-width reduction method
is elucidated. From a viewpoint of statistical mechanics, the BWR Ising model
has two-characteristic properties: an effective temperature and a slow
relaxation. These properties arise from the entropy effects of the auxiliary
spins, which are not present in the original Ising model.

  * •

To obtain the same dynamical process of the original Ising model, we propose
the setting parameters of the BWR Ising model for SA in which the temperature
schedule and inner loop are modified. The effectiveness of the proposed SA
parameters is evaluated using a dynamical process with square lattice random
Ising models. The dynamical process of the BWR Ising model is equivalent to
the original Ising model.

The rest of this paper is organized as follows. Section II introduces the bit-
width reduction method. Section III investigates the difference in dynamical
processes between the original Ising model, which has known properties, and
the BWR Ising model to clarify the dynamical properties of the BWR Ising
model. Section IV presents the statistical mechanics analysis results of the
BWR Ising model. Section V proposes BWR Ising model parameters for SA. The
experimental evaluations demonstrate that the dynamical processes of the BWR
Ising model and the original Ising model are almost the same. Sections VI and
VII demonstrate and discuss the numerical results, respectively. Section VIII
concludes with a summary of our study and future research directions. The
Appendices provide supplemental information for the derivation of the
statistical mechanics analysis for the BWR Ising model (Appendix A), the
effectiveness of the proposed SA parameter for large-size square lattice
systems (Appendix B) and various temperature schedules (Appendix C).

##  II method

###  II-A Bit-width reduction method

A previous study reported a bit-width reduction method, which added auxiliary
spins [32]. Herein the proposed bit-width reduction method modifies the
previous method for statistical mechanics analysis. Fig. 1 depicts the bit-
width reduction processes using the proposed method. Although the previous and
proposed methods have different coefficient assignments after bit-width
reduction, the substantive static properties are the same. Even with the
proposed method, the ground state of the BWR Ising model and the original
Ising model remain theoretically consistent.

Figure 1:  Examples of bit-width reduction using the proposed method. Arrows,
lines, circles, and squares represent the magnetic fields, interactions,
system spins, and auxiliary spins, respectively. (a) Bit-width reduction
process of the magnetic fields. Thick solid, dotted, and thin solid arrows
denote $h$, $h^{\prime}$, and $h^{\prime\prime}$, respectively. Lines between
system spin and auxiliary spins denote $|h^{\prime}|$. (b) Bit-width reduction
process of the interactions. Thick wavy, solid, dotted, and thin wavy lines
denote $J$, $J^{\prime}$, $|J^{\prime}|$, and $J^{\prime\prime}$,
respectively.

This subsection details the proposed method for bit-width reduction. The upper
and lower limits of the target $n$-bits coefficients are
$c^{\textrm{upper}}_{n}(=2^{n-1}-1)$ and
$c^{\textrm{lower}}_{n}(=-(2^{n-1}-1))$, respectively. First, we describe the
method to reduce the bit-width of the magnetic fields in the original Ising
model. We assume that the original Ising model includes a magnetic field $h$
acting on a spin $\sigma_{1}$ (Fig. 1, left). The spin consisting of the
original Ising model is called the “system spin” such as $\sigma_{1}$.
Applying the proposed method gives the BWR Ising model (Fig. 1, right). The
spins $s_{i}$ added by the proposed method are called “auxiliary spins.” The
bit-width of the magnetic field coefficient $h$ is reduced to $n$-bits as
follows:

Step 1:

Let $h^{\prime\prime}$ be the new magnetic field of $\sigma_{1}$, satisfying

$\displaystyle h^{\prime\prime}=\begin{cases}h-N_{\textrm{a}}\times
c^{\textrm{upper}}_{n},~{}0<h^{\prime\prime}\leqq
c^{\textrm{upper}}_{n}&(h>0)\\\ h-N_{\textrm{a}}\times
c^{\textrm{lower}}_{n},~{}c^{\textrm{lower}}_{n}\leqq
h^{\prime\prime}<0&(h<0)\end{cases}.$ (2)

Step 2:

Add $N_{\textrm{a}}$ auxiliary spins $s_{i}$ ($i=1,2,...,N_{\textrm{a}}$). Let
$h^{\prime}$ be the magnetic fields of all auxiliary spins $s_{i}$, where

$\displaystyle h^{\prime}=\begin{cases}c^{\textrm{upper}}_{n}&(h>0)\\\
c^{\textrm{lower}}_{n}&(h<0)\end{cases}.$ (3)

Step 3:

Introduce the interactions $|h^{\prime}|$ between $\sigma_{1}$ and all
auxiliary spins $s_{i}$.

Fig. 1 shows the scheme to reduce the bit-width of the interactions. We assume
that two spins $\sigma_{1}$ and $\sigma_{2}$ are connected by the interaction
$J$. The bit-width of the interaction coefficient $J$ is reduced to $n$-bits
as follows:

Step 1:

Let $J^{\prime\prime}$ be the new interaction between $\sigma_{1}$ and
$\sigma_{2}$, satisfying

$\displaystyle J^{\prime\prime}=\begin{cases}J-N_{\textrm{a}}\times
c^{\textrm{upper}}_{n},~{}0<J^{\prime\prime}\leqq
c^{\textrm{upper}}_{n}&(J>0)\\\ J-N_{\textrm{a}}\times
c^{\textrm{lower}}_{n},~{}c^{\textrm{lower}}_{n}\leqq
J^{\prime\prime}<0&(J<0)\end{cases}.$ (4)

Step 2:

Add $N_{\textrm{a}}$ auxiliary spins $s_{i}$ ($i=1,2,...,N_{\textrm{a}}$). Let
$J^{\prime}$ be the interactions between $\sigma_{1}$ all auxiliary spins
$s_{i}$, where

$\displaystyle J^{\prime}=\begin{cases}c^{\textrm{upper}}_{n}&(J>0)\\\
c^{\textrm{lower}}_{n}&(J<0)\end{cases}.$ (5)

Step 3:

Introduce the interactions $|J^{\prime}|$ between $\sigma_{2}$ and all
auxiliary spins $s_{i}$.

Figs. 2 and 2 show examples of the original Ising model and the BWR Ising
model obtained by applying the proposed method, respectively. In this case,
the bit-width of the coefficient is reduced from $4$-bits to $3$-bits. The
ground state of the original Ising model is ($\sigma_{1}$, $\sigma_{2}$,
$\sigma_{3}$)$=$($+1$, $+1$, $-1$). Similarly, the ground state of the BWR
Ising model is ($\sigma_{1}$, $\sigma_{2}$, $\sigma_{3}$, $s_{1}$, $s_{2}$,
$s_{3}$)$=$($+1$, $+1$, $-1$, $+1$, $-1$, $-1$). Focusing on the system spins,
the ground states of the original Ising model and the BWR Ising model are
clearly consistent.

Figure 2:  Example of bit-width reduction. Solid arrows, solid lines, circles,
and squares show the magnetic fields, their interactions, system spins, and
auxiliary spins, respectively. Orange and blue spins denote $+1$ and $-1$,
respectively. (a) Original Ising model ($4$-bits). (b) BWR Ising model after
applying the proposed method ($3$-bits).

###  II-B SA

SA is a meta-heuristic algorithm with a wide range of applications [48, 49,
50, 51]. During SA for $N$ spins Ising model given by (1), the following
procedures are performed:

Step 1:

Prepare a random initial spin state.

Step 2:

Set the initial temperature sufficiently high for the Hamiltonian.

Step 3:

Choose one spin from $N$ spins randomly.

Step 4:

Flip the chosen spin according to the transition probability $W(\Delta{E},T)$,
which depends on the temperature $T$ and the energy difference $\Delta E$.
Energy difference $\Delta E$ is defined by
$\Delta{E}=H_{\textrm{candidate}}-H_{\textrm{current}}$, where
$H_{\textrm{candidate}}$ is the energy of the candidate state in which the
chosen spin is flipped and $H_{\textrm{current}}$ is the energy of the current
state. Here, the transition probability, called the heat-bath method, is used
and is expressed as $W(\Delta{E},T)=[1+\exp(\Delta{E}/T)]^{-1}$.

Step 5:

Repeat Steps $3$–$4$ “inner loop” times. The inner loop is typically set to
the number of spins $N$, which is called one Monte Carlo Step (MCS).

Step 6:

Decrease the temperature $T$ and return to Step $3$.

Step 7:

Repeat Step $6$, “outer loop” times.

The Geman–Geman theorem guarantees that the ground state is ideally obtained
in SA when the temperature decreases sufficiently slow[52]. Notice that the
ground state may not be available and a lower-energy state (not the ground
state) may be obtained in a realistic time.

##  III Dynamical process with SA

In this study, the dynamical properties of the BWR Ising model were clarified
by comparing the dynamical process of the original Ising model to that of the
BWR Ising model under SA. We employed an Ising model on square $L\times L$
systems with periodic boundary conditions [53]. Here, we set $L=30$ and the
coefficients of magnetic fields and interactions to 7, that is, $h_{i}=7$ for
all $i$ and $J_{ij}=7$ for all nearest-neighbor pairs on square lattice $i,j$
in (1). The properties of the Ising model are well-known. In the ground state
of the original Ising model, all spins take $+1$, with an internal energy per
spin (i.e., energy density, $H/L^{2}$) of $-21$. In this demonstration, the
bit-width of the coefficient is reduced from $4$-bits to $3$\- or $2$-bits.

Table I shows the SA parameters. The initial temperature
$T_{\textrm{initial}}$ is set sufficiently high to permit the transition
between arbitrary states at the beginning of SA. The temperature schedule is
set to the power-law decay for every outer loop, which is given by
$T(t)=T_{\textrm{initial}}\times r^{t}$, where $r$ is the cooling rate and $t$
is the $t$-th outer loop. The outer loop and cooling rate $r$ is set to 100
($t=0-99$) and 0.97, respectively. Using these values, the final temperature
of SA becomes $2.451$, which is sufficiently low on the energy scale of the
original Ising model. The inner loop is set to the number of spins in the
Ising model ($1$ MCS).

TABLE I: SA parameters. Parameter | Value
---|---
Initial state | Random
Initial temperature ($T_{\textrm{initial}}$) | 50
Cooling rate ($r$) | 0.97
Outer loop | 100
Inner loop |  Number of spins (1 MCS)

Fig. 3 shows the experimental results. The energy density of both the original
and BWR Ising models was calculated using the number of system spins (i.e.,
$L^{2}$). The data were obtained from average and standard deviation of energy
density for ten simulations of SA. Although each BWR Ising model eventually
yields the ground state, the dynamical process significantly differs from that
of the original Ising model. Even at the steps of the outer loop, where the
energy density decreases in the original Ising model, it did not decrease in
the BWR Ising model. Similar results were obtained even for large-size square
lattice system ($L=40,50$) in Appendix B.

Figure 3: Dynamical processes of the original Ising model and the BWR Ising
model. The BWR Ising model ($2$-bits), the BWR Ising model ($3$-bits) and the
original Ising model ($4$-bits) are denoted by red squares, blue triangles and
black circles. Every plot is an average of ten runs. The error bars are
standard deviations.

##  IV Analysis of the bit-width reduced Ising model

To investigate the difference in the dynamical processes between the original
Ising model and the BWR Ising model, we analyzed the BWR Ising model from the
viewpoint of the microscopic mechanism: effective temperature and slow
relaxation.

###  IV-A Effective Temperature

Previous studies in statistical mechanics employed an Ising model with a
structure similar to the BWR Ising model (Fig. 1) [54, 55, 56]. It indicated
that the dynamical processes of the correlation function between the system
spins $\sigma_{1}$ and $\sigma_{2}$
$\lparen\langle\sigma_{1}\sigma_{2}\rangle\rparen$ of the temperature differ
from the Ising model with and without the auxiliary spins [54]. Here,
$\langle\cdot\rangle$ represents the expectation value. Therefore, we analyzed
the BWR Ising model by referencing the previous studies.

####  IV-A1 Magnetic fields

First, we considered the case where the bit-width of the magnetic fields is
reduced by adding $N_{\textrm{a}}$ auxiliary spins (Fig. 1). The effective
magnetic field $L_{\textrm{eff}}$ at temperature $T$ is defined as (see
Appendix A for a detailed derivation)

$\displaystyle\sum_{s_{i}=\pm
1}e^{-\beta{H}}=A(T)e^{L_{\mathrm{eff}}\sigma_{1}},$ (6)

where $\beta=1/T$ and $A(T)$ is an analytic function of $T$.

When $h>0$, the Hamiltonian of the BWR Ising model depicted on the right of
Fig. 1 is given by

$\displaystyle
H=-h^{\prime\prime}\sigma_{1}-h^{\prime}\left(\sum^{N_{\mathrm{a}}}_{i=1}s_{i}\right)(1+\sigma_{1}),$
(7)

and the effective magnetic field of the system spin is obtained as

$\displaystyle
L_{\mathrm{eff}}=\frac{h^{\prime\prime}}{T}+\frac{N_{\mathrm{a}}}{2}\log\left[\cosh\left(\frac{2h^{\prime}}{T}\right)\right].$
(8)

The effective temperature $T_{\textrm{eff}}=h/L_{\textrm{eff}}$ is given by

$\displaystyle
T_{\mathrm{eff}}=\frac{h}{\dfrac{h^{\prime\prime}}{T}+\dfrac{N_{\mathrm{a}}}{2}\log\left[\cosh\left(\dfrac{2h^{\prime}}{T}\right)\right]}.$
(9)

When $h<0$, the Hamiltonian depicted on the right of Fig. 1 is given by

$\displaystyle
H=-h^{\prime\prime}\sigma_{1}-h^{\prime}\left(\sum^{N_{\mathrm{a}}}_{i=1}s_{i}\right)(1-\sigma_{1}),$
(10)

and the effective magnetic field is obtained as

$\displaystyle
L_{\mathrm{eff}}=\frac{h^{\prime\prime}}{T}-\frac{N_{\mathrm{a}}}{2}\log\left[\cosh\left(\frac{2h^{\prime}}{T}\right)\right].$
(11)

The effective temperature is given by

$\displaystyle
T_{\mathrm{eff}}=\frac{h}{\dfrac{h^{\prime\prime}}{T}-\dfrac{N_{\mathrm{a}}}{2}\log\left[\cosh\left(\dfrac{2h^{\prime}}{T}\right)\right]}.$
(12)

####  IV-A2 Interactions

Next, we considered the case where the bit-width of the interactions is
reduced by adding $N_{\textrm{a}}$ auxiliary spins (Fig. 1).

The effective interaction $K_{\textrm{eff}}$ at temperature $T$ is defined as
(see Appendix A for a detailed derivation)

$\displaystyle\sum_{s_{i}=\pm
1}e^{-\beta{H}}=A(T)e^{K_{\mathrm{eff}}\sigma_{1}\sigma_{2}}.$ (13)

When $J>0$, the Hamiltonian of the BWR Ising model depicted on the right of
Fig. 1 is given by

$\displaystyle
H=-J^{\prime\prime}\sigma_{1}\sigma_{2}-J^{\prime}\left(\sum^{N_{\mathrm{a}}}_{i=1}s_{i}\right)(\sigma_{1}+\sigma_{2}),$
(14)

and the effective interaction between $\sigma_{1}$ and $\sigma_{2}$ is
obtained as

$\displaystyle
K_{\mathrm{eff}}=\frac{J^{\prime\prime}}{T}+\frac{N_{\mathrm{a}}}{2}\log\left[\cosh\left(\frac{2J^{\prime}}{T}\right)\right].$
(15)

The effective temperature $T_{\textrm{eff}}=h/K_{\textrm{eff}}$ is given by

$\displaystyle
T_{\mathrm{eff}}=\frac{J}{\dfrac{J^{\prime\prime}}{T}+\dfrac{N_{\mathrm{a}}}{2}\log\left[\cosh\left(\dfrac{2J^{\prime}}{T}\right)\right]}.$
(16)

When $J<0$, the Hamiltonian depicted on the right of Fig. 1 is given by

$\displaystyle
H=-J^{\prime\prime}\sigma_{1}\sigma_{2}-J^{\prime}\left(\sum^{N_{\mathrm{a}}}_{i=1}s_{i}\right)(\sigma_{1}-\sigma_{2}),$
(17)

and the effective interaction is obtained as

$\displaystyle
K_{\mathrm{eff}}=\frac{J^{\prime\prime}}{T}-\frac{N_{\mathrm{a}}}{2}\log\left[\cosh\left(\frac{2J^{\prime}}{T}\right)\right].$
(18)

The effective temperature is given by

$\displaystyle
T_{\mathrm{eff}}=\frac{J}{\dfrac{J^{\prime\prime}}{T}-\dfrac{N_{\mathrm{a}}}{2}\log\left[\cosh\left(\dfrac{2J^{\prime}}{T}\right)\right]}.$
(19)

Equations (9), (12), (16), and (19) indicate that the effective temperature
$T_{\textrm{eff}}$ differs from the temperature $T$ added to the Ising model.
Fig. 4 shows the effective temperature in the BWR Ising model, which was
determined by comparing the temperature $T$ used for SA in the previous
section and $T_{\textrm{eff}}$. Since $J=h=7$ is assumed, we set
$N_{\textrm{a}}=2$, $J^{\prime\prime}=1$, and $J^{\prime}=3$ for the
calculation to reduce the bit-width to $3$-bits. To reduce the bit-width of
coefficients to $2$-bits, we set $N_{\textrm{a}}=6$ and
$J^{\prime\prime}=J^{\prime}=1$. The temperature schedule of
$T_{\textrm{eff}}$ rapidly decreases at a temperature above that of $T$. This
suggests that the discrepancy between the temperature $T$ and
$T_{\textrm{eff}}$ affects the dynamical process.

Figure 4:  Effective temperature schedule of the original Ising model ($J$ or
$h=7$) and the BWR Ising model. Red squares, blue triangles, and black circles
denote the BWR Ising model ($2$-bits), BWR Ising model ($3$-bits), and
original Ising model ($4$-bits), respectively.

###  IV-B Characteristic time scale

The previous study reported that a slow relaxation occurs in the lattice of
frustrated systems with decorated spins [55]. This phenomenon is called
“entropic slowing down” and is due to the degrees of freedom distribution of
the decoration spins. The decorated lattice system has a similar structure to
the BWR Ising model when applying the proposed method. Therefore, we assumed
that an entropic slowing down appears in the BWR Ising model, and this
phenomenon influences the dynamical processes. Following [55], we determined
the number of states when the local configuration of the system spins is
fixed.

Figure 5:  (a) Parallel and (b) antiparallel system spin cases where the
original interactions are positive. Black, white, and gray symbols denote $+$,
$-$, and disordered spins, respectively. Circles and squares denote the system
and auxiliary spins, respectively.

First, we analyzed the case where the original interaction is positive (Fig.
5). Figs. 5 and 5 depict “the parallel state” (e.g.,
$(\sigma_{1},\sigma_{2})=(+,+)$) and “antiparallel state” of system spins
(e.g., $(\sigma_{1},\sigma_{2})=(+,-)$), respectively. Let $m$ denote the
number of auxiliary spins that represent the $+$ internal field when the
original interaction is positive. The energies of the parallel state and
antiparallel states are given by

$\displaystyle
E^{(N_{\mathrm{a}})}_{++}(m)=(2N_{\mathrm{a}}-4m)J^{\prime}-J^{\prime\prime},$
(20) $\displaystyle E^{(N_{\mathrm{a}})}_{+-}(m)=J^{\prime\prime}.$ (21)

We considered the probability distribution of the auxiliary spins at a
temperature $T$. In the parallel and antiparallel states, each probability of
$m$ up spins in the $N_{\textrm{a}}$ auxiliary spins are given by

$\displaystyle
Q^{(N_{\mathrm{a}})}_{++}(m)=\frac{\exp(-2\beta{J^{\prime}}N_{\mathrm{a}})}{(2\cosh{2\beta{J^{\prime}}})^{N_{\mathrm{a}}}}\binom{N_{\mathrm{a}}}{m}\exp(4\beta{J^{\prime}}m),$
(22) $\displaystyle
Q^{(N_{\mathrm{a}})}_{+-}(m)=\binom{N_{\mathrm{a}}}{m}\left(\frac{1}{2}\right)^{N_{\mathrm{a}}}.$
(23)

Next, we analyzed the case where the original interaction is negative. Figs. 6
and 6 depict “the parallel state” of system spins (e.g.,
$(\sigma_{1},\sigma_{2})=(-,-)$) and “antiparallel state” of system spins
(e.g., $(\sigma_{1},\sigma_{2})=(+,-)$), respectively. Let $n$ be the number
of auxiliary spins that represent the $+$ internal field when the original
interaction is negative. The energies of the parallel and antiparallel states
are given by

Figure 6:  (a) Parallel and (b) antiparallel system spin cases where the
original interactions are negative. Black, white, and gray symbols denote $+$,
$-$, and disordered spins, while the solid and dotted lines denote positive
and negative, respectively. Circles and squares denote the system spins and
the auxiliary spins, respectively.

$\displaystyle E^{(N_{\mathrm{a}})}_{--}(n)=-J^{\prime\prime},$ (24)
$\displaystyle
E^{(N_{\mathrm{a}})}_{+-}(n)=(2N_{\mathrm{a}}-4n)J^{\prime}+J^{\prime\prime}.$
(25)

Each probability of $n$ up spins in the $N_{\textrm{a}}$ auxiliary spins of
the parallel and antiparallel state are given by

$\displaystyle
R^{(N_{\mathrm{a}})}_{--}(n)=\binom{N_{\mathrm{a}}}{n}\left(\frac{1}{2}\right)^{N_{\mathrm{a}}},$
(26) $\displaystyle
R^{(N_{\mathrm{a}})}_{+-}(n)=\frac{\exp(-2\beta{J^{\prime}}N_{\mathrm{a}})}{(2\cosh{2\beta{J^{\prime}}})^{N_{\mathrm{a}}}}\binom{N_{\mathrm{a}}}{n}\exp(4\beta{J^{\prime}}n).$
(27)

Since $Q^{(N_{\textrm{a}})}_{+-}(m)$ and $R^{(N_{\textrm{a}})}_{--}(n)$ are
independent of the temperature, the equations are simple binomial
distributions and equivalent. In contrast, $Q^{(N_{\textrm{a}})}_{++}(m)$ and
$R^{(N_{\textrm{a}})}_{+-}(n)$ depend on temperature. They are maximized at
$N_{\textrm{a}}/2$ for high temperatures due to the entropy effect, whereas
they are maximized at nearly $N_{\textrm{a}}$ and $0$ for low temperatures.

Figure 7: Model for analyzing the flip probability of the central spin (free
spin) where the original interactions are (a) positive and (b) negative.
Black, white, and gray symbols denote $+$, $-$, and disordered spins,
respectively. Circles and squares denote the system and auxiliary spins,
respectively.

To estimate the transition probability of the system spins on the BWR Ising
model, we calculated the flip probability of the central spin shown in Fig. 7
following [55]. The central spin, which we refer to as “free spin,” is
surrounded by two up spins and two down spins. When $N_{\textrm{a}}=0$, the
flip probability of the free spin is $1/2$ in the Glauber dynamics [57].
However, when auxiliary spins are added by the proposed method
($N_{\textrm{a}}>0$), the flip probability becomes less than $1/2$ due to the
distribution of the surrounding auxiliary spins.

Fig. 7 shows the case where the original interaction is positive. The internal
field on the free spin is given by

$\displaystyle h(n_{1},n_{2})=2J^{\prime}(2N_{\mathrm{a}}-n_{1}-n_{2}),$ (28)

where $n_{1}$ and $n_{2}$ are the numbers of auxiliary spins representing the
$+$ internal field in the parallel and antiparallel state, respectively. In
the Glauber dynamics, the flip probability of free spin $P_{\textrm{flip}}$ is
given by

$P_{\mathrm{flip}}=\sum_{(n_{1},n_{2})}Q^{(2N_{\mathrm{a}})}_{++}(n_{1})Q^{(2N_{\mathrm{a}})}_{+-}(n_{2})\\\
\times\frac{1}{1+\exp[-2\beta{h}(n_{1},n_{2})]},$ (29)
$P_{\mathrm{flip}}=\frac{\exp(-4\beta{J^{\prime}}N_{\mathrm{a}})}{(4\cosh{2\beta{J^{\prime}}})^{2N_{\mathrm{a}}}}\sum_{(n_{1},n_{2})}\binom{2N_{\mathrm{a}}}{n_{1}}\binom{2N_{\mathrm{a}}}{n_{2}}\\\
\times\frac{\exp(4\beta{J^{\prime}}n_{1})}{1+\exp[-4\beta{J^{\prime}}(2N_{\mathrm{a}}-n_{1}-n_{2})]}.$
(30)

Similarly, in the case where the original interaction is negative (Fig. 7),
the internal field on a free spin is given by

$\displaystyle h(n_{3},n_{4})=-2J^{\prime}(n_{3}-n_{4}),$ (31)

where $n_{3}$ and $n_{4}$ are the numbers of auxiliary spins representing the
$+$ internal field in the antiparallel and parallel states, respectively. In
the Glauber dynamics, the $P_{\textrm{flip}}$ is given by

$P_{\mathrm{flip}}=\sum_{(n_{3},n_{4})}R^{(2N_{\mathrm{a}})}_{+-}(n_{3})R^{(2N_{\mathrm{a}})}_{--}(n_{4})\\\
\times\frac{1}{1+\exp[-2\beta{h}(n_{3},n_{4})]},$ (32)
$P_{\mathrm{flip}}=\frac{\exp(-4\beta{J^{\prime}}N_{\mathrm{a}})}{(4\cosh{2\beta{J^{\prime}}})^{2N_{\mathrm{a}}}}\sum_{(n_{3},n_{4})}\binom{2N_{\mathrm{a}}}{n_{3}}\binom{2N_{\mathrm{a}}}{n_{4}}\\\
\times\frac{\exp(4\beta{J^{\prime}}n_{3})}{1+\exp[4\beta{J^{\prime}}(n_{3}-n_{4})]}.$
(33)

The $P_{\textrm{flip}}$ is the same for arbitrary system spin combinations as
in (30) and (33) .

Figure 8: Flip probability of the free spin of the original Ising model
($J=7$) and the BWR Ising model for several temperatures. Red squares, blue
triangles, and black circles denote the BWR Ising model ($2$-bits), BWR Ising
model ($3$-bits), and original Ising model ($4$-bits), respectively.

Fig. 8 compares the flip probability of the free spin with the original and
BWR Ising models in the previous section. Although the probability is constant
with the number of auxiliary spins $N_{\textrm{a}}$ at high temperatures, it
changes significantly at low temperatures. Note that there is a limit to
$P_{\textrm{flip}}$ at low temperatures because the slow relaxation is caused
by the entropy effect. The value of the limit can be expressed as

$\displaystyle\lim_{T\to
0}P_{\mathrm{flip}}=\frac{1}{2}\left(\frac{1}{4}\right)^{N_{\mathrm{a}}}.$
(34)

A discrepancy in the flip probability occurs between the original and BWR
Ising models at low temperatures, which was not considered in the parameters
for SA in the previous section. This discrepancy likely affects the difference
in the dynamical processes. Note that the entropy effect does not occur in the
auxiliary spins for the magnetic fields.

##  V Proposed SA parameters

In the previous section, we analyzed the BWR Ising model using the proposed
method. The BWR Ising model has two characteristic properties: an effective
temperature and a slow relaxation. These properties are not present in the
original Ising model. In Section III, it was speculated that the dynamical
processes between the original and the BWR Ising model differ because the SA
is performed with the same SA parameters before and after bit-width reduction,
despite the variation in the statistical mechanics properties. This section
proposes SA parameters that consider the properties of the BWR Ising model and
evaluate the proposed SA parameters experimentally.

###  V-A How to modify the parameters

First, the temperature schedule is modified based on the effective temperature
$T_{\textrm{eff}}$ so that $T_{\textrm{eff}}$ is closer to the temperature $T$
of the original temperature schedule using (9), (12), (16), or (19). Fig. 9
shows the original temperature schedule used in Section III and the proposed
temperature schedules when the absolute value of the coefficient 0-7 is
reduced to $3$\- or $2$-bits.

Figure 9:  Proposed temperature schedules of the BWR Ising model for several
coefficients of the magnetic fields or interactions. Bit-widths are reduced to
(a) $3$-bits and (b) $2$-bits.

Next, the inner loop is modified based on the flip probability. To realize a
flip probability of the BWR Ising model closer to that of the original Ising
model, we define an effective relaxation time $\tau_{\textrm{eff}}$. According
to a previous study [55], $\tau_{\textrm{eff}}$ is given by

$\displaystyle\tau_{\mathrm{eff}}={P_{\mathrm{flip}}}^{-1}.$ (35)

Fig. 10 shows the relationship between temperature $T$ and
$\tau_{\textrm{eff}}$ of the original or the BWR Ising model.
$\tau_{\textrm{eff}}$ of the original Ising model is two from the definition
of the system shown in Fig. 7. The absolute value of coefficient 0-7 is
reduced to $3$\- or $2$-bits. Then $\tau_{\textrm{eff}}$ can be calculated by
(30) or (33), and (35).

Figure 10:  Effective time relaxation $\tau_{\textrm{eff}}$ of the BWR Ising
model for several coefficients of interactions. Bit-widths are reduced to (a)
$3$-bits width and (b) $2$-bits width.

The algorithm to modify the parameters is as follows:

Step 1:

Set the temperature schedule so that $T_{\textrm{eff}}$ is closer to the
original one.

Step 2:

Set the inner loop so that $\tau_{\textrm{eff}}$ of the BWR Ising model is
closer to that of the original Ising model. $\tau_{\textrm{eff}}$ is obtained
at each temperature determined in step 1. The proposed inner loop is set $1$
MCS$\times\tau_{\textrm{eff}}/2$ because the original inner loop corresponds
to $\tau_{\textrm{eff}}=2$.

###  V-B experimental evaluation

To evaluate the effectiveness of the proposed SA parameters, we experimentally
investigated the dynamical process of the energy density on the Ising model
used in Section III. We performed SA of the original Ising model with the SA
parameters described in Table I. For the BWR Ising model, the SA parameters
were changed from Table I to the proposed temperature schedule and the inner
loop explained in this section. We call this the “proposed SA parameters.” The
coefficients were reduced to $3$\- or $2$-bits. Fig. 11 shows the results. The
dynamical process of the BWR Ising model with the proposed SA parameters is
similar to that of the original Ising model. Similar effects were observed in
the large-size square lattice systems ($L=40,50$) in Appendix B.

Figure 11:  Dynamical process of the BWR Ising model with the proposed SA
parameters and the original Ising model with the original SA parameters. Red
squares, blue triangles, and black circles denote the BWR Ising model
($2$-bits), BWR Ising model ($3$-bits), and the original Ising model
($4$-bits), respectively. Every plot is the average of ten runs. The error
bars are standard deviations.

##  VI numerical results

We evaluated the applicability of the proposed SA parameters when the Ising
model has random coefficients. We compared the dynamical process of the
original and BWR Ising models using random Ising models [53]. We performed SA
of a square lattice system, where $L=30$. The coefficients of the magnetic
fields and interactions take integer values from $[-7,7]$ with equal
probabilities. Although $0$ was excluded for the interactions, it was included
for the magnetic fields. In these demonstrations, the bit-width of the
original Ising model was reduced from $4$-bits to $3$\- or $2$-bits.

Table I shows the SA parameters in this demonstration, except for the cooling
rate. The cooling rate $r$ was set such that the final temperature was equal
to $1$, i.e. $r=0.9612$. This condition ensures that the final temperature is
sufficiently small relative to the coefficients of the original Ising model.
See Appendix B for results using different types of temperature schedules.

In this demonstration, we performed SA with four types of SA parameters.

  * •

Original SA parameter
The unmodified SA parameters.

  * •

Modified temperature schedule (TS) SA parameter
Only the temperature schedule is modified based on the maximum absolute value
of the coefficients (i.e., $|J|,|h|=7$). A modified temperature schedule based
on the maximum absolute value shows the most gradual temperature decrease from
the lowest temperature (Fig. 9).

  * •

Modified inner loop (IL) SA parameter
Only the inner loop is modified based on the maximum absolute value of the
coefficients. Effective relaxation time $\tau_{\textrm{eff}}$ based on the
maximum absolute value shows the longest effective relaxation time (Fig. 10).

  * •

Proposed SA parameter
Both the temperature schedule and the inner loop are modified.

Fig. 12 shows the dynamical processes of the original and BWR Ising models for
each set of SA parameters. The data represent the average and standard
deviation of the energy density for ten SA simulations. By modifying the
temperature schedule, the energy density of the BWR Ising model becomes closer
to that of the original Ising model at the beginning of the iteration.
Applying an effective relaxation time to the inner loop prevents the energy
density of the BWR Ising model from terminating at a high value in the later
stages of the iteration. The dynamical processes of the BWR Ising model with
the proposed SA parameters are almost the same as that of the original Ising
model. However, a difference appears in the early stages of the iterations
when reducing to $3$-bits.

Figure 12:  Dynamical processes of the BWR Ising model and the original Ising
model. Black circles denote the original Ising model. Green triangles, blue
squares, orange diamonds, and red hexagons denote the BWR Ising model with the
original SA parameters, modified temperature schedule only, modified inner
loop only, and proposed SA parameters, respectively. Every plot is the average
of ten runs. The error bars are standard deviations.

Next, we performed SA with several sizes of Ising models [53] to evaluate the
problem size dependency of the proposed SA parameters. We set $L=10,20,30,40$,
or $50$ as the number of system spins. The coefficients and SA parameters are
as described above. Fig. 13 compares the energy densities for each Ising model
with a different size. Figs. 13 and 13 compare the energy densities at the end
of SA between the original and BWR Ising models using the original SA
parameters. All points are plotted above the diagonal, indicating that the BWR
Ising models have an inferior solution accuracy compared to the original Ising
model when using the original SA parameters. Figs. 13 and 13 compare the
energy densities at the end of SA between the original Ising model with the
original SA parameters and the BWR Ising model with the proposed SA
parameters. All points are plotted below or on the diagonal, indicating that
the BWR Ising model is the same or superior to the original Ising model when
using the proposed SA parameters. These results imply that this feature is
independent of the system size, at least for the range considered in this
study. Additionally, it was demonstrated that the dynamical process is also
independent of the system size (Appendix B). These results suggest that our
proposed SA parameters exhibit robustness to spin size in the square lattice
system.

Figure 13:  Comparisons of the energy densities for several sizes of Ising
models between the BWR Ising model and the original Ising model with the
original SA parameters when the bit-widths are reduced to (a) $3$-bits or (c)
$2$-bits, or the BWR Ising model with the proposed SA parameters and the
original Ising model with the original SA parameters when the bit-widths are
reduced to (b) $3$-bits or (d) $2$-bits. Every plot is the average of ten
runs. The error bars are standard deviations.

##  VII Discussion

The dynamical process of the BWR Ising model with the proposed SA parameters
is almost the same as that of the original Ising model, although the early
stages of the iterations differ when reduced to $3$-bits (Fig. 12). The
difference is attributed to the coefficient used as the basis for modifying
the SA parameters. In the previous section, we used the maximum absolute
values of the original Ising model as the basis. However, there is a gap
between the temperature schedule based on the maximum absolute value of the
coefficient and the other coefficients (Fig. 9). The temperature schedule
based on the maximum absolute value of the coefficient has the lowest
temperature in the early stages of the outer loop iterations and the slowest
temperature decrease.

Additionally, the coefficients of the Ising model used in the previous section
were generated uniformly at random. In many scenarios, the coefficients are
not seven. That is, the proportion of coefficients (i.e., $|J|$, $|h|=6$ to
$4$) with a large gap from the temperature schedule based on the maximum
absolute value of coefficient when reduced to $3$-bits is relatively high
(Fig. 9). This indicates that the proposed temperature schedule is set
excessively low and slow for many coefficients. Consequently, the dynamical
process of the BWR Ising model in the early stages of outer loop iterations
differs from that of the original Ising model.

The solution accuracy of the BWR Ising models with the proposed SA parameters
is slightly superior to that of the original Ising model (Figs. 13 and 13).
This is attributed to the coefficient used as the basis for modifying the SA
parameters. The effective relaxation time $\tau_{\textrm{eff}}$ based on the
maximum absolute value of the coefficient is larger than that based on the
other coefficients at a low temperature (Fig. 10). Because the proposed inner
loop is set by an excess $\tau_{\textrm{eff}}$ based on an excessively low and
slow temperature schedule, for many coefficients, the probability of reaching
thermal equilibrium at each temperature is higher and the solution accuracy is
improved.

##  VIII Conclusion and future work

The dynamical process with SA is compared between the original Ising model and
the BWR Ising model by applying the proposed method using square lattice
systems. Because the dynamical process of the BWR Ising model significantly
differs from that of the original Ising model, we analyzed the BWR Ising model
from the viewpoint of statistical mechanics. The BWR Ising model with the
addition of auxiliary spins has two-characteristic properties not present in
the original Ising model: an effective temperature and a slow relaxation.
Considering the analytical results, we proposed SA parameters for the BWR
Ising model. Our results demonstrate that the dynamical processes of the BWR
Ising model with the proposed SA parameters are close to that of the original
Ising model.

We expect that the bit-width reduction method and our parameter modification
method will effectively solve the Ising model on an implemented Ising machine
with a bit-width limitation of coefficients. However, these methods are not
efficient in terms of computation time. The computation time increases as the
number of auxiliary spins increases. The number of auxiliary spins
$N_{\mathrm{a}}$ to be added by the proposed method per one coefficient is
represented by:

$\displaystyle N_{\mathrm{a}}\simeq\begin{cases}2^{\lparen
n_{\mathrm{original}}-n_{\mathrm{BWR}}\rparen}&(n_{\mathrm{BWR}}>2)\\\
2^{\lparen n_{\mathrm{original}}-1\rparen}-2&(n_{\mathrm{BWR}}=2),\end{cases}$
(36)

where $n_{\textrm{original}}$ and $n_{\textrm{BWR}}$ are the bit-widths of the
original Ising model and the BWR Ising model, respectively. The increased
computation time is directly related to the inner loop. The modified inner
loop is set to $1$ MCS$\times\tau_{\textrm{eff}}$. The MCS increases with the
number of auxiliary spins because MCS is the total number of spins of the
Ising model. $\tau_{\textrm{eff}}$ increases with the number of auxiliary
spins and temperature. The modified temperature schedule is lower than that of
the original temperature schedule due to the increased number of auxiliary
spins. Therefore, the computation time increases.

A method has been proposed that combines the shift method with the bit-width
reduction method using auxiliary spins, to reduce the number of auxiliary
spins required [58]. However, this approach leads to a different ground state
than the original Ising model. One alternative method to maintain the ground
state while mitigating the increase in computation time is to expand the bit-
width that can be input to the Ising machine. However, even if the bit-width
of the Ising machine cannot be increased due to hardware limitations, the
number of the MCS can be reduced if the auxiliary spins can be flipped
simultaneously similar to CMOS annealing [59]. This approach should be
considered in the future. Additionally, we plan to investigate other
implemented algorithms of Ising machines such as quantum annealing [56].

##  Appendix A Analysis of the entropic effects caused by auxiliary spin

This appendix provides detailed derivations of the definitions for the
effective magnetic field $L_{\textrm{eff}}$ and the effective interaction
$K_{\textrm{eff}}$ (see (6) in the main text).

####  A-1 Magnetic fields

First, consider the case of a bit-width reduction for a magnetic field. The
Ising model is assumed to be the same as that in Section IV-A1 (Fig. 1). To
derive the definition of the effective magnetic field, the expectation value
of the system spin $\langle\sigma_{1}\rangle$ of the original Ising model and
that of the BWR Ising model are matched. The expectation value of the original
Ising model $\langle\sigma_{1}\rangle_{\textrm{original}}$ is given by

$\displaystyle\langle\sigma_{1}\rangle_{\mathrm{original}}=\sum_{\sigma_{1}=\pm
1}\sigma_{1}\times P=\frac{\sum_{\sigma_{1}=\pm
1}\sigma_{1}e^{L{\sigma_{1}}}}{\sum_{\sigma_{1}=\pm 1}e^{L{\sigma_{1}}}},$
(37)

where $L=h/T$ and $P$ is the probability distribution at temperature $T$. $P$
is given by

$\displaystyle P=\frac{e^{-\beta H}}{\sum_{\sigma_{i}=\pm 1}e^{-\beta H}}.$
(38)

The expectation value of the BWR Ising
model$\langle\sigma_{1}\rangle_{\textrm{BWR}}$ is similarly given by

$\displaystyle\langle\sigma_{1}\rangle_{\mathrm{BWR}}=\frac{\sum_{\sigma_{1}=\pm
1}\sigma_{1}e^{L_{\mathrm{eff}}{\sigma_{1}}}}{\sum_{\sigma_{1}=\pm
1}e^{L_{\mathrm{eff}}{\sigma_{1}}}},$ (39)

where $L_{\textrm{eff}}=h/T_{\textrm{eff}}$ and $T_{\textrm{eff}}$ is the
effective temperature. To match (37) and (39), the auxiliary spins in (39) are
partially summed to obtain the marginal probability for $\sigma_{1}$.
Therefore, the effective magnetic field $L_{\textrm{eff}}$ is defined by
tracing the auxiliary spins $s_{i}$ to give (6) in the main text.

####  A-2 Interactions

In the case of a bit-width reduction of the interactions, the Ising model
shown in Fig. 1 is assumed. The expectation values of the original Ising model
$\langle\sigma_{1}\sigma_{2}\rangle_{\textrm{original}}$ and the BWR Ising
model $\langle\sigma_{1}\sigma_{2}\rangle_{\textrm{BWR}}$ are given by

$\displaystyle\langle\sigma_{1}\sigma_{2}\rangle_{\mathrm{original}}=\frac{\sum_{\sigma_{i}=\pm
1}\sigma_{1}\sigma_{2}e^{K{\sigma_{1}\sigma_{2}}}}{\sum_{\sigma_{i}=\pm
1}e^{K{\sigma_{1}\sigma_{2}}}},$ (40)
$\displaystyle\langle\sigma_{1}\sigma_{2}\rangle_{\mathrm{BWR}}=\frac{\sum_{\sigma_{i}=\pm
1}\sigma_{1}\sigma_{2}e^{K_{\mathrm{eff}}{\sigma_{1}\sigma_{2}}}}{\sum_{\sigma_{i}=\pm
1}e^{K_{\mathrm{eff}}{\sigma_{1}\sigma_{2}}}},$ (41)

where $K=J/T$ and $K_{\textrm{eff}}=J/T_{\textrm{eff}}$, respectively. To
match (40) and (41), the effective interaction $K_{\textrm{eff}}$ is defined
by tracing the auxiliary spins $s_{i}$, which yields (13) in the main text.

##  Appendix B Dynamical process for large-size square lattice system

In this appendix, to evaluate the problem size dependency of the dynamical
process, we performed SA on a large-size square lattice system ($L=40,50$)
[53] with several types of SA parameters.

We first investigated the dynamical process of the square lattice system with
all coefficients set to 7, as described in Section III. The dynamical
processes are shown in Fig. 14. The results revealed that even in large-size
square lattice systems, the dynamical processes of the original Ising model
and the BWR Ising model differ when using the original SA parameters, shown in
Table I. However, by using the proposed SA parameters described in Section V,
the dynamical processes of the original Ising model and the BWR Ising model
became almost the same.

Figure 14:  Dynamical processes of the original Ising model with the original
SA parameters and the BWR Ising model with (a), (c) the original SA
parameters, or (b), (d) the proposed SA parameters. Square lattice system size
are (a), (b) $L=40$ and (c), (d) $L=50$. The BWR Ising model ($2$-bits), the
BWR Ising model ($3$-bits) and the original Ising model ($4$-bits) are denoted
by red squares, blue triangles and black circles. Every plot is an average of
ten runs. The error bars are standard deviations.

Next, we investigated the dynamical process of a large-size square lattice
system with randomly assigned coefficients, as described in Section VI. We
also used the SA parameters as described in the same Section VI. Figure 15
shows the dynamical processes. As mentioned in Section VI, the dynamical
processes of the BWR Ising model with the proposed SA parameters became closer
to that of the original Ising model.

Figure 15:  Dynamical processes of the original Ising model and the BWR Ising
model reduced to (a), (c) $3$-bits or (b), (d) $2$-bits. Square lattice system
sizes are (a), (b) $L=40$ and (c), (d) $L=50$. Black circles denote the
original Ising model. Green triangles, blue squares, orange diamonds, and red
hexagons denote the BWR Ising model with the original SA parameters, modified
temperature schedule only, modified inner loop only, and proposed SA
parameters, respectively. Every plot is the average of ten runs. The error
bars are standard deviations.

##  Appendix C Performance of the proposed method for different temperature
schedules

This appendix evaluates the performance of the proposed SA parameters when SA
is performed at various temperature schedules. Table II and Fig. 16 show the
temperature schedules. The initial state, outer loop, and inner loop were set
to random, $100$ ($n=0-99$), and 1 MCS, respectively, as described in Sections
III and VI. These parameters were adjusted so that the final temperature was
around one.

Figure 16: Different types of temperature schedules. The $n$ in the temperature schedule represents the $n$-th outer loop. TABLE II: Different temperature schedules in the SA. The $n$ in the temperature schedule represents the $n$-th outer loop. Temperature schedule |  Initial Temperature ($T_{\textrm{initial}}$) |  Final Temperature
---|---|---
$T(n)=T_{\textrm{initial}}\times 0.9612^{n}$ | 50 | 0.994
$T(n)=T_{\textrm{initial}}\times 0.394n$ | 40 | 0.994
$T(n)=T_{\textrm{initial}}/n$ | 50 | 0.990
$T(n)=T_{\textrm{initial}}/\sqrt{n}$ | 50 | 1.005
$T(n)=T_{\textrm{initial}}/\ln{(1+n)}$ | 40 | 1.506

We performed SA for a square lattice system with $L=20$. The coefficients of
the magnetic fields and interactions take values under the same conditions as
those described in Section VI.

Figure 17:  Dynamical processes of the BWR Ising model and the original Ising
model. Black circles denote the original Ising model. The $n$ in the
temperature schedule represents the $n$-th outer loop. Green triangles and
orange squares denote the BWR Ising model with the original SA parameters when
the bit-widths are reduced to a $3$-bit width or $2$-bit width, respectively.
Blue triangles and red squares denote the BWR Ising model with the proposed SA
parameters when the bit-widths are reduced to $3$-bits or $2$-bits,
respectively. Every plot is the average of ten runs. The error bars are
standard deviations.

Fig. 17 shows the dynamical processes of the BWR Ising model with the proposed
or original parameters and the original Ising model with the original
parameters using a different temperature schedule. For all temperature
schedules, the dynamical processes and the energy densities at the end of SA
of the BWR Ising model with the original SA parameters differed from that of
the original Ising model. In contrast, the BWR Ising model with the proposed
SA parameters showed a similar performance as the original Ising model.

## Acknowledgment

This article is based on the results obtained from a project, JPNP16007,
commissioned by the New Energy and Industrial Technology Development
Organization (NEDO). The computation in this work has been partially done
using the facilities of the Supercomputer Center, the Institute for Solid
State Physics, the University of Tokyo. S. T. was supported in part JSPS
KAKENHI (Grant Numbers JP21K03391, JP23H05447) and JST Grant Number
JPMJPF2221. Human Biology-Microbiome-Quantum Research Center (Bio2Q) is
supported by World Premier International Research Center Initiative (WPI),
MEXT, Japan.

## References

  * [1] A. Lucas, “Ising formulations of many NP problems,” _Front. Phys._ , vol. 2, no. 5, pp. 1–15, 2014.
  * [2] S. Tanaka, R. Tamura, and B. K. Chakrabarti, _Quantum Spin Glasses, Annealing and Computation_. Cambridge University Press, 2017.
  * [3] M. W. Johnson, M. H. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, “Quantum annealing with manufactured spins,” _Nature_ , vol. 473, no. 7346, pp. 194–198, 2011.
  * [4] A. Askarzadeh, L. dos Santos Coelho, C. E. Klein, and V. C. Mariani, “A population-based simulated annealing algorithm for global optimization,” in _2016 IEEE international conference on systems, man, and cybernetics (SMC)_. IEEE, 2016, pp. 004 626–004 633.
  * [5] C. Yamaguchi, “Proposal of a checking parameter in the simulated annealing method applied to the spin glass model,” _Comput. Phys. Commun._ , vol. 199, pp. 47–52, 2016.
  * [6] 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 J. Solid-State Circuits_ , vol. 51, no. 1, pp. 303–309, 2016.
  * [7] M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and H. G. Katzgraber, “Physics-inspired optimization for quadratic unconstrained problems using a digital annealer,” _Front. Phys._ , vol. 7, no. 48, pp. 1–14, 2019.
  * [8] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, O. Tadanaga, H. Takenouchi, K. Aihara, K.-i. Kawarabayashi, K. Inoue, S. Utsunomiya, and H. Takesue, “A coherent Ising machine for 2000-node optimization problems,” _Science_ , vol. 354, no. 6312, pp. 603–606, 2016.
  * [9] H. Goto, K. Tatsumura, and A. R. Dixon, “Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems,” _Sci. Adv._ , vol. 5, no. 4, pp. 1–9, 2019.
  * [10] M. Maezawa, G. Fujii, M. Hidaka, K. Imafuku, K. Kikuchi, H. Koike, K. Makise, S. Nagasawa, H. Nakagawa, M. Ukibe, and S. Kawabata, “Toward practical-scale quantum annealing machine for prime factoring,” _J. Phys. Soc. Jpn._ , vol. 88, no. 6, p. 061012, 2019.
  * [11] 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 J. Solid-State Circuits_ , vol. 56, no. 1, pp. 165–178, 2020.
  * [12] N. Mohseni, P. L. McMahon, and T. Byrnes, “Ising machines as hardware solvers of combinatorial optimization problems,” _Nat. Rev. Phys._ , vol. 4, no. 6, pp. 363–379, 2022.
  * [13] H. Neven, V. S. Denchev, M. Drew-Brook, J. Zhang, W. G. Macready, and G. Rose, “NIPS 2009 demonstration: Binary classification using hardware implementation of quantum annealing,” _Quantum_ , vol. 4, 2009.
  * [14] M. H. Amin, “Searching for quantum speedup in quasistatic quantum annealers,” _Phys. Rev. A_ , vol. 92, p. 052323, 2015.
  * [15] M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and R. Melko, “Quantum Boltzmann machine,” _Phys. Rev. X_ , vol. 8, p. 021050, 2018.
  * [16] D. O’Malley, V. V. Vesselinov, B. S. Alexandrov, and L. B. Alexandrov, “Nonnegative/binary matrix factorization with a D-Wave quantum annealer,” _PloS one_ , vol. 13, no. 12, p. e0206653, 2018.
  * [17] K. Kitai, J. Guo, S. Ju, S. Tanaka, K. Tsuda, J. Shiomi, and R. Tamura, “Designing metamaterials with quantum annealing and factorization machines,” _Phys. Rev. Res._ , vol. 2, p. 013319, Mar 2020.
  * [18] T. Inoue, Y. Seki, S. Tanaka, N. Togawa, K. Ishizaki, and S. Noda, “Towards optimization of photonic-crystal surface-emitting lasers via quantum annealing,” _Opt. Express_ , vol. 30, no. 24, pp. 43 503–43 512, 2022\.
  * [19] K. Endo, Y. Matsuda, S. Tanaka, and M. Muramatsu, “A phase-field model by an Ising machine and its application to the phase-separation structure of a diblock polymer,” _Sci. Rep._ , vol. 12, no. 1, p. 10794, 2022.
  * [20] G. Rosenberg, P. Haghnegahdar, P. Goddard, P. Carr, K. Wu, and M. L. De Prado, “Solving the optimal trading trajectory problem using a quantum annealer,” _IEEE J. Sel. Top. Signal Processing_ , vol. 10, no. 6, pp. 1053–1060, 2016\.
  * [21] K. Tanahashi, S. Takayanagi, T. Motohashi, and S. Tanaka, “Application of Ising machines and a software development for Ising machines,” _J. Phys. Soc. Jpn._ , vol. 88, no. 6, p. 061010, 2019.
  * [22] A. Perdomo-Ortiz, N. Dickson, M. Drew-Brook, G. Rose, and A. Aspuru-Guzik, “Finding low-energy conformations of lattice protein models by quantum annealing,” _Sci. Rep._ , vol. 2, no. 571, 2012.
  * [23] F. Neukart, G. Compostella, C. Seidel, D. von Dollen, S. Yarkoni, and B. Parney, “Traffic flow optimization using a quantum annealer,” _Front. ICT_ , vol. 4, p. 29, 2017.
  * [24] H. Irie, G. Wongpaisarnsin, M. Terabe, A. Miki, and S. Taguchi, “Quantum annealing of vehicle routing problem with time, state and capacity,” in _International Workshop on Quantum Technology and Optimization Problems_. Springer, 2019, pp. 145–156.
  * [25] S. Bao, M. Tawada, S. Tanaka, and N. Togawa, “An approach to the vehicle routing problem with balanced pick-up using Ising machines,” in _2021 International Symposium on VLSI Design, Automation and Test (VLSI-DAT)_ , 2021, pp. 1–4.
  * [26] ——, “Multi-day travel planning using Ising machines for real-world applications,” in _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_ , 2021, pp. 3704–3709.
  * [27] Y. Mukasa, T. Wakaizumi, S. Tanaka, and N. Togawa, “An Ising machine-based solver for visiting-route recommendation problems in amusement parks,” _IEICE Trans. Info. Syst._ , vol. 104, no. 10, pp. 1592–1600, 2021.
  * [28] S. Naito, Y. Hasegawa, Y. Matsuda, and S. Tanaka, “ISAAQ: Ising machine assisted quantum compiler,” _arXiv preprint arXiv:2303.02830_ , 2023.
  * [29] S. Izawa, K. Kitai, S. Tanaka, R. Tamura, and K. Tsuda, “Continuous black-box optimization with an Ising machine and random subspace coding,” _Phys. Rev. Res._ , vol. 4, no. 2, p. 023062, 2022.
  * [30] Y. Seki, R. Tamura, and S. Tanaka, “Black-box optimization for integer-variable problems using Ising machines and factorization machines,” _arXiv preprint arXiv:2209.01016_ , 2022.
  * [31] S. Tanaka, Y. Matsuda, and N. Togawa, “Theory of Ising machines and a common software platform for Ising machines,” in _2020 25th Asia and South Pacific Design Automation Conference (ASP-DAC)_ , 2020, pp. 659–666.
  * [32] D. Oku, M. Tawada, S. Tanaka, and N. Togawa, “How to reduce the bit-width of an Ising model by adding auxiliary spins,” _IEEE Trans. Comput._ , vol. 71, no. 1, pp. 223–234, 2022.
  * [33] M. Kowalsky, T. Albash, I. Hen, and D. A. Lidar, “3-regular three-xorsat planted solutions benchmark of classical and quantum heuristic optimizers,” _Quantum Sci. Technol._ , vol. 7, no. 2, p. 025008, 2022.
  * [34] M. Booth, S. P. Reinhardt, and A. Roy, “Partitioning optimization problems for hybrid classical/quantum execution,” in _D-Wave Technical Report Series 14-1006A-A_ , 2017, https://docs.ocean.dwavesys.com/projects/qbsolv/en/latest/index.html.
  * [35] H. Karimi and G. Rosenberg, “Boosting quantum annealer performance via sample persistence,” _Quantum Inf. Process_ , vol. 16, no. 7, pp. 1–24, 2017.
  * [36] H. Karimi, G. Rosenberg, and H. G. Katzgraber, “Effective optimization using sample persistence: A case study on quantum annealers and various monte carlo optimization methods,” _Phys. Rev. E_ , vol. 96, no. 4, p. 043312, 2017.
  * [37] S. Okada, M. Ohzeki, M. Terabe, and S. Taguchi, “Improving solutions by embedding larger subproblems in a d-wave quantum annealer,” _Sci. Rep._ , vol. 9, no. 2098, pp. 1–10, 2019.
  * [38] H. Irie, H. Liang, T. Doi, S. Gongyo, and T. Hatsuda, “Hybrid quantum annealing via molecular dynamics,” _Sci. Rep._ , vol. 11, no. 1, p. 8426, 2021.
  * [39] Y. Atobe, M. Tawada, and N. Togawa, “Hybrid annealing method based on subqubo model extraction with multiple solution instances,” _IEEE Trans. Comput._ , vol. 71, no. 10, pp. 2606–2619, 2022.
  * [40] S. Kikuchi, N. Togawa, and S. Tanaka, “Hybrid optimization method using simulated-annealing-based Ising machine and quantum annealer,” _arXiv preprint arXiv:2306.06589_ , 2023.
  * [41] V. Choi, “Minor-embedding in adiabatic quantum computation: I. the parameter setting problem,” _Quantum Inf. Process._ , vol. 7, no. 5, pp. 193–209, 2008\.
  * [42] ——, “Minor-embedding in adiabatic quantum computation: II. minor-universal graph design,” _Quantum Inf. Process._ , vol. 10, no. 3, pp. 343–353, 2011.
  * [43] J. Cai, W. G. Macready, and A. Roy, “A practical heuristic for finding graph minors,” _arXiv preprint arXiv:1406.2741_ , 2014.
  * [44] T. Boothby, A. D. King, and A. Roy, “Fast clique minor generation in chimera qubit connectivity graphs,” _Quantum Inf. Process._ , vol. 15, no. 1, pp. 495–508, 2016.
  * [45] A. Zaribafiyan, D. J. Marchand, and S. S. Changiz Rezaei, “Systematic and deterministic graph minor embedding for cartesian products of graphs,” _Quantum Inf. Process._ , vol. 16, no. 5, pp. 1–26, 2017.
  * [46] D. Oku, K. Terada, M. Hayashi, M. Yanaoka, S. Tanaka, and N. Togawa, “A fully-connected Ising model embedding method and its evaluation for CMOS annealing machines,” _IEICE Trans. Info. Syst._ , vol. E102.D, no. 9, pp. 1696–1706, 2019.
  * [47] T. Shirai, S. Tanaka, and N. Togawa, “Guiding principle for minor-embedding in simulated-annealing-based Ising machines,” _IEEE Access_ , vol. 8, pp. 210 490–210 502, 2020.
  * [48] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” _Science_ , vol. 220, no. 4598, pp. 671–680, 1983.
  * [49] D. S. Johnson, C. R. Aragon, L. A. McGeoch, and C. Schevon, “Optimization by simulated annealing: An experimental evaluation; part I, graph partitioning,” _Oper. Res._ , vol. 37, no. 6, pp. 865–892, 1989.
  * [50] ——, “Optimization by simulated annealing: an experimental evaluation; part II, graph coloring and number partitioning,” _Oper. Res._ , vol. 39, no. 3, pp. 378–406, 1991.
  * [51] S. V. Isakov, I. N. Zintchenko, T. F. Rønnow, and M. Troyer, “Optimised simulated annealing for Ising spin glasses,” _Comput. Phys. Commun._ , vol. 192, pp. 265–271, 2015.
  * [52] S. Geman and D. Geman, “Stochastic relaxation, gibbs distributions, and the bayesian restoration of images,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. PAMI-6, no. 6, pp. 721–741, 1984.
  * [53] “Dataset of square lattice systems,” [Online]. Available: https://github.com/Shu-Tanaka-Group/bit-width-reduced-Ising-model, accessed on: July, 2023.
  * [54] S. Miyashita, S. Tanaka, and M. Hirano, “Nonmonotonic relaxation in systems with reentrant-type interaction,” _J. Phys. Soc. Jpn._ , vol. 76, no. 8, p. 083001, 2007.
  * [55] S. Tanaka and S. Miyashita, “Mechanism of slow relaxation due to screening effect in a frustrated system,” _J. Phys. Soc. Jpn._ , vol. 78, no. 8, p. 084002, 2009.
  * [56] S. Tanaka and S. Miyashita, “Nonmonotonic dynamics in a frustrated Ising model with time-dependent transverse field,” _Phys. Rev. E_ , vol. 81, no. 5, p. 051138, 2010.
  * [57] R. J. Glauber, “Time-dependent statistics of the Ising model,” _J. Math. Phys._ , vol. 4, no. 2, pp. 294–307, 1963.
  * [58] Y. Yachi, M. Tawada, and N. Togawa, “An efficient combined bit-width reducing method for Ising models,” _IEICE Trans. Inf. Syst._ , vol. 106, no. 4, pp. 495–508, 2023.
  * [59] T. Okuyama, M. Hayashi, and M. Yamaoka, “An Ising computer based on simulated quantum annealing by path integral Monte Carlo method,” _2017 IEEE International Conference on Rebooting Computing, ICRC 2017 - Proceedings_ , vol. 2017-Janua, pp. 1–6, 2017.

|  Shuta Kikuchi received the B.Eng. and M.Eng. degrees from the Waseda
University in 2017 and 2019, respectively. He is currently pursuing a Ph.D.
degree in applied physics at Keio University. His research interests include
Ising machine, statistical mechanics, and quantum annealing. He is a member of
the JPS.
---|---
|  Nozomu Togawa (Member, IEEE) received the B.Eng., M.Eng., and Dr.Eng.
degrees in electrical engineering from the Waseda University, Tokyo, Japan, in
1992, 1994, and 1997, respectively. He is currently a Professor with the
Department of Computer Science and Communications Engineering, Waseda
University. His research interests include quantum computation and integrated
system design. He is a member of ACM, IEICE, and IPSJ.
---|---
|  Shu Tanaka received the B.Sci. degree from the Tokyo Institute of
Technology, in 2003, and the M.Sci. and Dr.Sci. degrees from The University of
Tokyo, in 2005 and 2008, respectively. He is currently an Associate Professor
with the Department of Applied Physics and Physico-Informatics, Keio
University and a Core Director with Human Biology-Microbiome-Quantum Research
Center (Bio2Q), Keio University. His research interests include quantum
annealing, Ising machine, statistical mechanics, and materials science. He is
a member of the JPS.
---|---


# Effect of solenoid lens field on electron beam emittance

Kazuaki Togawa<EMAIL_ADDRESS>RIKEN SPring-8 Center, 1-1-1 Kouto, Sayo,
Hyogo 679-5148, Japan

###### Abstract

In an injector system of an X-ray free electron laser (XFEL), solenoid lenses
are typically used to confine low-emittance electron beams to low-energy
region below a few MeV. Because non-thermionic emittance at such a low-energy
region is easily deteriorated by nonlinear electromagnetic fields, it is
important to determine the properties of a solenoid lens on electron beam
emittance in the design of XFEL injectors. We derived an approximate solution
to emittance growth due to lens aberration by a paraxial approximation. It was
found that the derivative of the longitudinal magnetic field strongly affects
beam emittance, and its growth is proportional to the fourth power of the beam
radius. Various properties of the beam can be analyzed as long as the
longitudinal magnetic field distribution is prepared using a simulation or
measurement. In this study, a theoretical procedure to obtain the emittance
growth in the solenoid lens is introduced and the design considerations of the
solenoid lens of the SACLA injector are described.

###### pacs:

07.55.Db, 41.60.Cr

##  I Introduction

X-ray free electron lasers (XFELs) require an extremely fine electron beam
whose emittance is below 1 $\mu$m to lase high-intensity coherent X-ray beams
[1]. Two types of electron guns are used currently for XFELs: photo-cathode
radio-frequency and thermionic-cathode high-voltage guns [2, 3]. Because the
electron beams must be accelerated and transported in the low-energy region
from 0 to a few MeV in either case, the non-thermionic emittances are easily
deteriorated by nonlinear electromagnetic fields such as space charge and
accelerating fields.

At a low-energy region, a series of solenoid coils or lenses are frequently
used to confine the electron beam because of the strong space charge force
that results in the beam spreading out instantly. It is common that the
solenoid lens has an intrinsic aberration property. This property can
determine the limit of imaging resolution in electron microscopes [4]. In
electron accelerators, it causes emittance growth that deteriorates XFEL
output power. Therefore, it is important to evaluate the basic properties of a
solenoid lens to minimize emittance growth in electron injector design. In
this study, an approximate solution of emittance growth in the solenoid lens
is derived and the design consideration for the low-emittance injector of the
Japanese XFEL facility, SACLA [5], is described.

##  II Emittance growth in solenoid lens

###  II.1 Equation of motion in axisymmetric magnetic field

The equation of motion of an electron in an axisymmetric magnetic field with
cylindrical coordinates is given by

$\displaystyle\gamma m_{e}\frac{dv_{r}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\gamma m_{e}v_{\theta}^{2}}{r}-ev_{\theta}B_{z},$ (1)
$\displaystyle\gamma m_{e}\frac{d(rv_{\theta})}{dt}$ $\displaystyle=$
$\displaystyle-er(v_{z}B_{r}-v_{r}B_{z}),$ (2)

where $\gamma$ is the relativistic mass factor, $m_{e}$ is the electron mass,
$e$ is the elementary charge, $v_{\theta}$ and $v_{z}$ are the rotational and
longitudinal components of velocity, respectively, and $B_{r}$ and $B_{z}$ are
the radial and longitudinal components of the magnetic field, respectively.
Here, the rotational component $B_{\theta}$ that is induced by the beam
current is ignored and $\gamma$ is assumed to be constant.

The velocity increase in the radial direction is obtained by integrating
Eq.(1):

$\displaystyle\Delta
v_{r}=\int\biggl{(}\frac{v_{\theta}^{2}}{r}-\frac{e}{\gamma
m_{e}}v_{\theta}B_{z}\biggr{)}dt.$ (3)

Using $dz=v_{z}dt$ and $v_{z}\approx const.$, Eq.(3) can be rewritten as

$\displaystyle\Delta
v_{r}=\frac{1}{v_{z}}\int\biggl{(}\frac{v_{\theta}^{2}}{r}-\frac{e}{\gamma
m_{e}}v_{\theta}B_{z}\biggr{)}dz.$ (4)

Next, we derive the relation between the velocity in the rotational direction
$v_{\theta}$ and rotational component of vector potential $A_{\theta}$. From
$\overrightarrow{\rm B}=rot\overrightarrow{\rm A}$, $B_{r}$ and $B_{z}$ can be
expressed as

$\displaystyle B_{r}$ $\displaystyle=$ $\displaystyle-\frac{\partial
A_{\theta}}{\partial z},$ (5) $\displaystyle B_{z}$ $\displaystyle=$
$\displaystyle\frac{1}{r}\frac{\partial(rA_{\theta})}{\partial r},$ (6)

where $\partial A_{z}/\partial\theta=\partial A_{r}/\partial\theta=0$ is used
because the magnetic field is axisymmetrical. Using Eq.(5) and (6), Eq.(2)
becomes

$\displaystyle\gamma
m_{e}\frac{d(rv_{\theta})}{dt}=e\frac{d(rA_{\theta})}{dt}.$ (7)

When $A_{\theta}$ is zero at the initial point, i.e., the magnetic field does
not penetrate a beam emitter cathode, $v_{\theta}$ can be expressed as

$\displaystyle v_{\theta}=\frac{e}{\gamma m_{e}}A_{\theta}.$

Substituting $v_{\theta}$ in Eq.(4), $\Delta v_{r}$ becomes

$\displaystyle\Delta v_{r}=\frac{1}{v_{z}}\biggl{(}\frac{e}{\gamma
m_{e}}\biggr{)}^{2}\int\biggl{(}\frac{A_{\theta}^{2}}{r}-A_{\theta}B_{z}\biggr{)}dz.$

As a result, the variation of the electron diverging angle $\Delta r^{\prime}$
after passing through the solenoid field becomes:

$\displaystyle\Delta r^{\prime}=\frac{\Delta
v_{r}}{v_{z}}=\biggl{(}\frac{e}{\beta\gamma
m_{e}c}\biggr{)}^{2}\int\biggl{(}\frac{A_{\theta}^{2}}{r}-A_{\theta}B_{z}\biggr{)}dz,$
(8)

where $\beta=v_{z}/c$.

###  II.2 Paraxial expansion of magnetic field

To obtain the paraxial approximation of $A_{\theta}$ and $B_{z}$, we start
with the Maxwell equation:

$\displaystyle\nabla^{2}\overrightarrow{\rm A}=0.$

The rotational component can be written explicitly as

$\displaystyle\frac{1}{r}\frac{\partial}{\partial r}\biggl{(}r\frac{\partial
A_{\theta}}{\partial r}\biggr{)}+\frac{\partial^{2}A_{\theta}}{\partial
z^{2}}-\frac{A_{\theta}}{r^{2}}=0,$ (9)

where $\partial A_{\theta}/\partial\theta=0$ is used because the magnetic
field is axisymmetrical. Moreover, $A_{\theta}$ can be expanded as

$\displaystyle A_{\theta}(r,z)=\sum^{\infty}_{n=0}a_{n}(z)r^{n}.$ (10)

Substituting Eq.(10) in Eq.(9), we obtain the relation of the coefficients
$a_{n}$:

$\displaystyle\sum^{\infty}_{n=0}(1-n^{2})a_{n}r^{n-2}=\sum^{\infty}_{n=0}a_{n}^{(2)}r^{n},$
(11)

where $a_{n}^{(k)}$ shows $k$th derivative of the coefficient:

$\displaystyle a_{n}^{(k)}=\frac{d^{k}a_{n}}{dz^{k}}.$

By comparing each coefficients of $r^{n}$ in Eq.(11), we obtain

$\displaystyle a_{0}=0,$ $\displaystyle[1-(n+2)^{2}]a_{n+2}=a_{n}^{(2)}.$ (12)

Eq.(12) can be rewritten as

$\displaystyle a_{n+2}=-\frac{a_{n}^{(2)}}{(n+1)(n+3)}.$

From $a_{0}=0$, all of the coefficients with even index numbers become zero.
The coefficients with odd index numbers are

$\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle a_{1},$ $\displaystyle
a_{3}$ $\displaystyle=$ $\displaystyle-\frac{a_{1}^{(2)}}{2\cdot 4},$
$\displaystyle a_{5}$ $\displaystyle=$
$\displaystyle-\frac{a_{3}^{(2)}}{4\cdot 6}=(-1)^{2}\frac{1}{2\cdot 4\cdot
4\cdot 6}a_{1}^{(4)},$ $\displaystyle a_{7}$ $\displaystyle=$
$\displaystyle-\frac{a_{5}^{(2)}}{6\cdot 8}=(-1)^{3}\frac{1}{2\cdot 4\cdot
4\cdot 6\cdot 6\cdot 8}a_{1}^{(6)},$ $\displaystyle\cdot$ $\displaystyle\cdot$
$\displaystyle\cdot$ $\displaystyle a_{n}$ $\displaystyle=$
$\displaystyle\frac{(-1)^{\frac{n-1}{2}}}{2^{n-1}\left[\left(\frac{n-1}{2}\right)!\right]^{2}\frac{n+1}{2}}a_{1}^{(n-1)}\
\ \ \ (n=3,5,7,\cdot\cdot\cdot).$

Because $n$ is an odd number, it can be replaced by $2n+1$:

$\displaystyle a_{2n+1}=\frac{(-1)^{n}}{2^{2n}n!(n+1)!}a_{1}^{(2n)}\ \ \ \
(n=1,2,3,\cdot\cdot\cdot).$

Consequently, the vector potential $A_{\theta}$ is expressed using the
following polynomial:

$\displaystyle
A_{\theta}(r,z)=\sum^{\infty}_{n=0}\frac{(-1)^{n}}{2^{2n}n!(n+1)!}a_{1}^{(2n)}(z)r^{2n+1}.$
(13)

The longitudinal magnetic field $B_{z}$ can be expressed as a polynomial using
Eq.(6) and Eq.(13):

$\displaystyle
B_{z}(r,z)=\sum^{\infty}_{n=0}\frac{(-1)^{n}}{2^{2n}(n!)^{2}}2a_{1}^{(2n)}(z)r^{2n}.$

Here, $a_{1}$ gives the longitudinal component of the z-axis magnetic field:

$\displaystyle B_{z}(0,z)$ $\displaystyle=$ $\displaystyle
B_{z0}(z)=2a_{1}(z).$

Using $B_{z0}$, we obtain the final expressions of $A_{\theta}$ and $B_{z}$:

$\displaystyle A_{\theta}(r,z)$ $\displaystyle=$
$\displaystyle\sum^{\infty}_{n=0}\frac{(-1)^{n}}{n!(n+1)!}B_{z0}^{(2n)}(z)\left(\frac{r}{2}\right)^{2n+1},$
(14) $\displaystyle B_{z}(r,z)$ $\displaystyle=$
$\displaystyle\sum^{\infty}_{n=0}\frac{(-1)^{n}}{(n!)^{2}}B_{z0}^{(2n)}(z)\left(\frac{r}{2}\right)^{n}.$
(15)

###  II.3 Emittance growth owing to nonlinear focusing

To evaluate beam focusing properties of the solenoid lens, $A_{\theta}$ and
$B_{z}$ are approximated using the first and second terms in Eq.(14) and
Eq.(15), respectively:

$\displaystyle A_{\theta}(r,z)$ $\displaystyle=$
$\displaystyle\frac{1}{2}rB_{z0}(z)-\frac{1}{16}r^{3}B^{\prime\prime}_{z0}(z),$
(16) $\displaystyle B_{z}(r,z)$ $\displaystyle=$ $\displaystyle
B_{z0}(z)-\frac{1}{4}r^{2}B^{\prime\prime}_{z0}(z),$ (17)

where $B^{\prime\prime}_{z0}$ denotes the second derivative $B_{z0}^{(2)}$.
Using Eq.(16) and Eq.(17), the integrand in Eq.(8) becomes

$\displaystyle-\frac{r}{4}B_{z0}^{2}(z)+\frac{r^{3}}{8}B_{z0}(z)B_{z0}^{\prime\prime}(z)-\frac{3r^{5}}{256}B_{z0}^{\prime\prime
2}(z).$

The last term with the 5th order of $r$ is negligible. Substituting them in
Eq.(8), the increase in the beam diverging angle $\Delta r^{\prime}$ can be
expressed as

$\displaystyle\Delta r^{\prime}=\left(\frac{e}{\beta\gamma
m_{e}c}\right)^{2}\left(-\int\frac{B_{z0}^{2}(z)}{4}dz\cdot
r+\int\frac{B_{z0}(z)B_{z0}^{\prime\prime}(z)}{8}dz\cdot r^{3}\right).$ (18)

Using the following relation

$\displaystyle\int B_{z0}^{\prime 2}(z)dz+\int
B_{z0}(z)B_{z0}^{\prime\prime}(z)dz=[B_{z0}(z)B_{z0}^{\prime}(z)]_{-\infty}^{\infty}=0,$

Eq.(18) becomes

$\displaystyle\Delta r^{\prime}=-\frac{1}{8}\left(\frac{e}{\beta\gamma
m_{e}c}\right)^{2}\left(2\int B_{z0}^{2}(z)dz\cdot r+\int B_{z0}^{\prime
2}(z)dz\cdot r^{3}\right).$ (19)

In an ordinary case, a field region created by the thin lens is short and the
a radial position of electron $r$ does not change significantly. Therefore, we
assume that $r$ is constant; however $r^{\prime}$ changes in the solenoid
field. Linear beam focusing is obtained from the first term, whereas the
second term causes nonlinear beam focusing, i.e., lens aberration. It should
be noted that the sign of the nonlinear term is always similar to that of the
linear one. This means that the lens aberration cannot be compensated using
another lens.

The focal length of the solenoid lens $f$ can be estimated using the first
term of Eq.(19):

$\displaystyle\frac{1}{f}=-\frac{\Delta
r^{\prime}}{r}=\frac{1}{4}\left(\frac{e}{\beta\gamma m_{e}c}\right)^{2}\int
B_{z0}^{2}(z)dz.$ (20)

In a low energy case, using $(\beta\gamma)^{2}\approx 2eV/(m_{e}c^{2})$, the
focal length is approximated as

$\displaystyle\frac{1}{f}=\frac{e}{8m_{e}V}\int B_{z0}^{2}(z)dz,$

where $V$ is the beam voltage. This is a familiar formula commonly described
in electron optics textbooks [6].

Next, we consider the emittance growth owing to nonlinear focusing. The
normalized rms emittance increase in the horizontal (or vertical) direction is
given by

$\displaystyle\Delta\varepsilon_{n,rms}^{x}=\frac{\beta\gamma}{2}\sqrt{<r^{2}><\Delta
r^{\prime 2}>-<r\cdot\Delta r^{\prime}>^{2}},$ (21)

where $<>$ denotes the mean values weighted by the current density. Each mean
value can be calculated using Eq.(19):

$\displaystyle<r^{2}>$ $\displaystyle=$ $\displaystyle\frac{\int
J(r)r^{2}dS}{\int J(r)dS}=\frac{\int r^{3}dr}{\int rdr}=\frac{r^{2}}{2},$ (22)
$\displaystyle<\Delta r^{\prime 2}>$ $\displaystyle=$ $\displaystyle\frac{\int
J(r)\Delta r^{\prime 2}dS}{\int
J(r)dS}=\frac{\int(2k_{1}r+k_{2}r^{3})^{2}rdr}{\int
rdr}=2k_{1}^{2}r^{2}+\frac{4}{3}k_{1}k_{2}r^{4}+\frac{1}{4}k_{2}^{2}r^{6},$
(23) $\displaystyle<r\cdot\Delta r^{\prime}>$ $\displaystyle=$
$\displaystyle\frac{\int J(r)r\Delta r^{\prime}dS}{\int
J(r)dS}=\frac{\int(2k_{1}r+k_{2}r^{3})r^{2}dr}{\int
rdr}=k_{1}r^{2}+\frac{1}{3}k_{2}r^{4},$ (24)

where

$\displaystyle k_{1}$ $\displaystyle=$
$\displaystyle-\frac{1}{8}\left(\frac{e}{\beta\gamma m_{e}c}\right)^{2}\int
B_{z0}^{2}(z)dz,$ $\displaystyle k_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{8}\left(\frac{e}{\beta\gamma m_{e}c}\right)^{2}\int
B_{z0}^{\prime 2}(z)dz,$

and the beam current density is assumed to be uniform ($J(r)=const.$).
Substituting Eq.(22)-(24) in Eq.(21), the emittance increase becomes

$\displaystyle\Delta\varepsilon_{n,rms}^{lens}=\frac{\beta\gamma}{2}\frac{|k_{2}|}{6\sqrt{2}}r^{4}.$

Finally, we obtain the approximative expression of the emittance growth in the
solenoid lens:

$\displaystyle\Delta\varepsilon_{n,rms}^{lens}=\frac{1}{96\sqrt{2}\beta\gamma}\left(\frac{e}{m_{e}c}\right)^{2}\int\left(\frac{dB_{z0}(z)}{dz}\right)^{2}dz\cdot
r^{4}.$ (25)

This equation suggests the following two important features: 1) Emittance
increases intensely when beam size becomes large. Beam size should be
sufficiently small to maintain low emittance. In addition, the law that
emittance growth is proportional to the fourth power of the beam radius
appears in an electric field [7]. 2) The derivative of the longitudinal
magnetic field affects beam emittance. Therefore, the magnetic field should
vary gradually along the z-axis. This can be achieved by selecting a large
aperture of the solenoid yoke.

##  III Solenoid lens design considerations for a low-emittance injector

In the SACLA injector system and SCSS test accelerator [8], which was upgraded
to SCSS+ [9], a total of ten solenoid lenses are used to focus and transport
the beam whose energy range is 0.5-1 MeV. Therefore, it is essential to
minimize the non-thermionic emittance growth owing to the solenoid fields. To
determine the shape of solenoid yoke, three types of solenoid lenses were
analyzed. Type (a), which was chosen for the SCSS test accelerator and SACLA,
has a large pole piece aperture of yoke. Type (b) has a small aperture, which
was preferable in creating a strong magnetic field in a narrow space. Type (c)
was composed of two solenoid lenses whose magnetic field polarity were
opposite. The opposite fields can cancel the rotation around the beam axis.
The cross section of the solenoid lenses are shown in Fig. 1.

Figure 1:  Cross section of the solenoid lenses with: (a) a large pole piece
aperture of yoke, (b) a small pole piece aperture of yoke, and (c) a double
lens whose magnetic field polarity are opposite. The dimensions are described
in mm. The SCSS test accelerator and SACLA adopted type (a) with outer yoke
diameters of 190 mm and 304 mm, respectively.

Fig. 2 shows the longitudinal magnetic field distribution on the z-axis
calculated using the POISSON code [10]. The coil currents for each type were
assumed to be 2297 A$\cdot$turn for (a), 1554 A$\cdot$turn for (b), and 4236
A$\cdot$turn for (c). A focul length of 0.3 m for a 500 keV beam was obtained
from these currents, which was determined using Eq.(20). When the magnetic
flux penetrates the cathode, normalized rms emittance increases by

$\displaystyle\Delta\varepsilon_{n,rms}^{mag}=\frac{eB_{zc}r_{c}^{2}}{8m_{e}c},$

where $B_{zc}$ is the longitudinal magnetic field on the cathode and $r_{c}$
is the cathode radius. The large pole piece aperture causes a wide magnetic
field distribution. Therefore, the distance between the cathode and first
solenoid lens must be long sufficiently to certify the small field. In SACLA,
the distance was selected 135 mm.

Figure 2: Longitudinal magnetic field on z-axis calculated using the POISSON
code.

The $r$-$\Delta r^{\prime}$ phase space profiles after passing through the
solenoid lens can be calculated using Eq.(19). The integral of the square of
the magnetic field and derivative were calculated using field data shown in
Fig. 2. It is assumed that the beam energy is 500 keV and the radial position
of electron in the field is constant. It is clearly shown in Fig. 3 that the
nonlinearity of the profile curve was distinctive for small aperture and
double lenses.

Figure 3: Phase space profiles of electron beams after passing through the
solenoid lenses.

The emittance growth was calculated using Eq.(25) and is shown in Fig. 4. It
can be seen that the emittance growth is significantly suppressed by enlarging
the pole piece aperture. In the SACLA injector, the beam diameter is
controlled to less than $\sim\phi$5 mm. In this case, the emittance growth is
less than 0.01 $\mu$m per each lens. Therefore, it is predicted that the total
emittance growth should not exceed 0.2 $\mu$m.

Figure 4: Emittance growths owing to solenoid lens aberration as a function of
the beam radius.

The abovementioned emittance growth is caused by geometric distortion of the
linearly distributed phase space profile, i.e., a non-thermal phenomenon. A
beam extracted from a hot cathode has a thermal emittance
$\varepsilon_{n,rms}^{th}$ determined by its temperature and radius. The net
emittance is given by their root-sum-square of:

$\displaystyle\varepsilon_{n,rms}^{net}=\sqrt{{\varepsilon_{n,rms}^{th}}^{2}+{\Delta\varepsilon_{n,rms}^{lens}}^{2}}.$

The CeB6 gun generates a beam with a low thermal emittance of 0.6 $\mu$m [3].
Therefore, the increment of net emittance is $\sim$0.03 $\mu$m. This value is
negligibly small for the XFEL injector. It should be noted that
$\varepsilon_{n,rms}^{net}$ represents the sliced beam emittance, which is
important in achieving FEL gain saturation for XFEL. Increment of the
projected emittance owing to time dependent energy dispersion in the bunched
beam is not mentioned here.

According to the preceding discussion, a smaller beam size is favorable in
minimizing lens aberration; however, the space charge force in the beam
becomes stronger. The equation of motion of an electron in a space charge
dominated beam is described by

$\displaystyle\frac{d^{2}r}{dz^{2}}=\frac{2I(r)}{I_{0}(\beta\gamma)^{3}}\frac{1}{r},$
$\displaystyle I(r)=\int_{0}^{r}J(r_{1})2{\pi}r_{1}dr_{1},$

where $I_{0}$ is Alfén current (17 kA) and $J(r)$ is a beam current density at
a radial position $r$. If the current density is nonuniform, the nonlinearity
of the beam diverging angle $r^{\prime}$ is reinforced by the strong space
charge. In particular, it becomes severe after the sub-harmonic buncher, where
the beam is bunched up. Therefore, the beam size must be tuned as a compromise
between minimizing the lens aberration and avoiding emittance growth owing to
the nonlinear space charge.

##  IV Summary

In this study, we derived the approximative solution of an emittance growth
owing to the solenoid lens aberration. It was found that the emittance
develops strongly when the incident beam has a large diameter. In addition, it
was evident that the derivative of longitudinal magnetic field affected the
emittance. Using derived formulae, the solenoid lens for the SACLA injector
was designed to sustain low beam emittance. Therefore, it can be concluded
that emittance growth owing to lens aberration is negligibly small for SACLA.

###### Acknowledgements.

The author would like to thank Dr. T. Shintake of Okinawa Institute of Science
and Technology Graduate University, Dr. A. D. Yeremian of SLAC National
Laboratory, and V. A. Goryashko of Uppsala University for valuable discussion
during this research.

## References

  * McNeil and Thompson [2010] B. W. J. McNeil and N. R. Thompson, Nat. Photonics 4, 814 (2010).
  * Akre et al. [2008] R. Akre, D. Dowell, P. Emma, J. Frisch, S. Gilevich, G. Hays, P. Hering, R. Iverson, C. Limborg-Deprey, H. Loos, et al., Phys. Rev. ST Accel. Beams 11, 030703 (2008), URL https://link.aps.org/doi/10.1103/PhysRevSTAB.11.030703.
  * Togawa et al. [2007] K. Togawa, T. Shintake, T. Inagaki, K. Onoe, T. Tanaka, H. Baba, and H. Matsumoto, Phys. Rev. ST Accel. Beams 10, 020703 (2007), URL https://link.aps.org/doi/10.1103/PhysRevSTAB.10.020703.
  * Tonomura [1987] A. Tonomura, Rev. Mod. Phys. 59, 639 (1987), URL https://link.aps.org/doi/10.1103/RevModPhys.59.639.
  * Ishikawa et al. [2012] T. Ishikawa, H. Aoyagi, T. Asaka, Y. Asano, N. Azumi, T. Bizen, H. Ego, K. Fukami, T. Fukui, Y. Furukawa, et al., Nat. Photon. 6, 540 (2012).
  * Pierce [1954] J. R. Pierce, _Theory and Design of Electron Beams_ (D. Van Nostrand Company, Inc., Princeton, New Jersey, 1954), 2nd ed.
  * Opanasenko et al. [2017] A. Opanasenko, V. Mytrochenko, V. Zhaunerchyk, and V. A. Goryashko, Phys. Rev. Accel. Beams 20, 053401 (2017), URL https://link.aps.org/doi/10.1103/PhysRevAccelBeams.20.053401.
  * Shintake et al. [2008] T. Shintake, H. Tanaka, T. Hara, T. Tanaka, K. Togawa, M. Yabashi, Y. Otake, T. Inagaki, S. Inoue, M. Ishii, et al., Nat. Photon. 2, 555 (2008).
  * Owada et al. [2018] S. Owada, K. Togawa, T. Inagaki, T. Hara, T. Tanaka, Y. Joti, T. Koyama, K. Nakajima, H. Ohashi, Y. Senba, et al., J. Synchrotron Rad. 25, 282 (2018).
  * Menzel and Stokes [1987] M. T. Menzel and H. K. Stokes, LA-UR 87$-$115, Los Alamos National Laboratory (1987).