Dynamical process of a bit-width reduced Ising model with simulated annealing

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.


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 op-timization 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 where σ i is the spin on the vertex i ∈ V and has a value of (+1) or (−1). h i is the magnetic field on the vertex i ∈ V , and J ij is the interaction on the edge (i, j) ∈ E. In this paper, we assume that interactions and magnetic fields are integer constants. VOLUME  Approaches such as meta-heuristics and Ising machines have been developed to solve combinatorial optimization problems [3]- [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]- [16], material design [17]- [19], portfolio optimization [20], [21], protein folding [22], traffic optimization [23]- [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]- [40] and the topology [41]- [47], few studies have been devoted to overcoming the bit-width limitation.

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 A. BIT-WIDTH REDUCTION METHOD
A previous study reported a bit-width reduction method, which added auxiliary spins [32]. Herein the proposed bitwidth 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.
This subsection details the proposed method for bit-width reduction. The upper and lower limits of the target n-bits coefficients are c upper n (= 2 n−1 −1) and c lower n (= −(2 n−1 −1)), respectively. First, we describe the method to reduce the bitwidth 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 σ 1 (Fig. 1(a), left). The spin consisting of the original Ising model is called the ''system spin'' such as σ 1 . Applying the proposed method gives the BWR Ising model ( Fig. 1(a), right). The spins s i added by the proposed 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 ′ , and h ′′ , respectively. Lines between system spin and auxiliary spins denote |h ′ |. (b) Bit-width reduction process of the interactions. Thick wavy, solid, dotted, and thin wavy lines denote J, J ′ , |J ′ |, and J ′′ , respectively. 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 ′′ be the new magnetic field of σ 1 , satisfying Step 2:Add N a auxiliary spins s i (i = 1, 2, ..., N a ). Let h ′ be the magnetic fields of all auxiliary spins s i , where Step 3:Introduce the interactions |h ′ | between σ 1 and all auxiliary spins s i . Fig. 1(b) shows the scheme to reduce the bit-width of the interactions. We assume that two spins σ 1 and σ 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 ′′ be the new interaction between σ 1 and σ 2 , satisfying Step 2:Add N a auxiliary spins s i (i = 1, 2, ..., N a ). Let J ′ be the interactions between σ 1 all auxiliary spins s i , where Step 3:Introduce the interactions |J ′ | between σ 2 and all auxiliary spins s i . Figs. 2(a) and 2(b) 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 (σ 1 , σ 2 , σ 3 )=(+1, +1,  SA is a meta-heuristic algorithm with a wide range of applications [48]- [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 (∆E, T ), which depends on the temperature T and the energy difference ∆E. Energy difference ∆E is defined by ∆E = H candidate − H current , where H candidate is the energy of the candidate state in which the chosen spin is flipped and H current is the energy of the current state. Here, the transition probability, called the heat-bath method, is used and is expressed as W (∆E, T ) = [1 + exp(∆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 × L systems with periodic boundary conditions [53]. Here, we set L = 30 and the coefficients of magnetic fields and interactions to 7, VOLUME 11, 2023 In this demonstration, the bit-width of the coefficient is reduced from 4-bits to 3or 2bits. Table 1 shows the SA parameters. The initial temperature T 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 initial × 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). 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.

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.

A. EFFECTIVE TEMPERATURE
Previous studies in statistical mechanics employed an Ising model with a structure similar to the BWR Ising model ( Fig. 1(b)) [54]- [56]. It indicated that the dynamical processes of the correlation function between the system spins σ 1 and σ 2 (⟨σ 1 σ 2 ⟩) of the temperature differ from the Ising model with and without the auxiliary spins [54]. Here, ⟨·⟩ represents the expectation value. Therefore, we analyzed the BWR Ising model by referencing the previous studies.

1) Magnetic fields
First, we considered the case where the bit-width of the magnetic fields is reduced by adding N a auxiliary spins ( Fig. 1(a)). The effective magnetic field L eff at temperature T is defined as (see Appendix A for a detailed derivation) where β = 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

(a) is given by
and the effective magnetic field of the system spin is obtained as The effective temperature T eff = h/L eff is given by When h < 0, the Hamiltonian depicted on the right of Fig. 1

(a) is given by
and the effective magnetic field is obtained as The effective temperature is given by 2) Interactions Next, we considered the case where the bit-width of the interactions is reduced by adding N a auxiliary spins ( Fig. 1(b)). The effective interaction K eff at temperature T is defined as (see Appendix A for a detailed derivation) When J > 0, the Hamiltonian of the BWR Ising model depicted on the right of Fig. 1

(b) is given by
and the effective interaction between σ 1 and σ 2 is obtained as The effective temperature T eff = h/K eff is given by When J < 0, the Hamiltonian depicted on the right of Fig. 1

(b) is given by
and the effective interaction is obtained as The effective temperature is given by Equations (9), (12), (16), and (19) indicate that the effective temperature T 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 eff . Since J = h = 7 is assumed, we set N a = 2, J ′′ = 1, and J ′ = 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 a = 6 and J ′′ = J ′ = 1. The temperature schedule of T eff rapidly decreases at a temperature above that of T . This suggests that  the discrepancy between the temperature T and T eff affects the dynamical process.

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.
First, we analyzed the case where the original interaction is positive (Fig. 5). Figs. 5(a) and 5(b) depict ''the parallel state'' (e.g., (σ 1 , σ 2 ) = (+, +)) and ''antiparallel state'' of system spins (e.g., (σ 1 , σ 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 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 a auxiliary spins are given by Next, we analyzed the case where the original interaction is negative. Figs. 6(a) and 6(b) depict ''the parallel state'' of system spins (e.g., (σ 1 , σ 2 ) = (−, −)) and ''antiparallel state'' of system spins (e.g., (σ 1 , σ 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 Each probability of n up spins in the N a auxiliary spins of the parallel and antiparallel state are given by Since Q 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 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 a > 0), the flip probability becomes less than 1/2 due to the distribution of the surrounding auxiliary spins. Fig. 7(a) shows the case where the original interaction is positive. The internal field on the free spin is given by 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 flip is given by .
Similarly, in the case where the original interaction is negative (Fig. 7(b)), the internal field on a free spin is given by 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 flip is given by .
The P flip is the same for arbitrary system spin combinations as in (30) and (33) . 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 a at high temperatures, it changes significantly at low temperatures. Note that there is a limit to P flip at low temperatures because the slow relaxation is caused by the entropy effect. The value of the limit can be expressed as 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.

A. HOW TO MODIFY THE PARAMETERS
First, the temperature schedule is modified based on the effective temperature T eff so that T 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 3or 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 τ eff . According to a previous study [55], τ eff is given by Fig. 10 shows the relationship between temperature T and τ eff of the original or the BWR Ising model. τ 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 3or 2-bits. Then τ eff can be calculated by (30) or (33), and (35).
The algorithm to modify the parameters is as follows: Step 1:Set the temperature schedule so that T eff is closer to the original one.
Step 2:Set the inner loop so that τ eff of the BWR Ising model is closer to that of the original Ising model. τ eff is obtained at each temperature determined in step 1. The proposed inner loop is set 1 MCS×τ eff /2 because the original inner loop corresponds to τ eff = 2.

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 1. For the BWR Ising model, the SA parameters were changed from Table 1 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 3or 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.

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 3or 2-bits. Table 1 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 τ 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.
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(a) and 13(c) 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(b) and 13(d) 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.

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 3bits is relatively high (Fig. 9(a)). 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(b) and 13(d)). This is attributed to the coefficient used as the basis for modifying the SA parameters. The effective relaxation time τ 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 τ 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 a to be added by the proposed method per one coefficient is represented by: where n original and n 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×τ eff . The MCS increases with the number of auxiliary spins because MCS is the total number of spins of the Ising model. τ 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 eff and the effective interaction K eff (see (6) in the main text).

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(a)). To derive the definition of the effective magnetic field, the expectation value of the system spin ⟨σ 1 ⟩ of the original Ising model and that of the BWR Ising model are matched. The expectation value of the original Ising model ⟨σ 1 ⟩ original is given by where L = h/T and P is the probability distribution at temperature T . P is given by The expectation value of the BWR Ising model⟨σ 1 ⟩ BWR is similarly given by where L eff = h/T eff and T eff is the effective temperature. To match (37) and (39), the auxiliary spins in (39) are partially summed to obtain the marginal probability for σ 1 . Therefore, the effective magnetic field L eff is defined by tracing the auxiliary spins s i to give (6) in the main text.

2) Interactions
In the case of a bit-width reduction of the interactions, the Ising model shown in Fig. 1(b) is assumed. The expectation values of the original Ising model ⟨σ 1 σ 2 ⟩ original and the BWR Ising model ⟨σ 1 σ 2 ⟩ BWR are given by where K = J /T and K eff = J /T eff , respectively. To match (40) and (41), the effective interaction K 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 1. 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.
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.

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 2 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. 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. 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 Indus-