1 Introduction↩︎

The slowing of Moore’s law in terms of advancement in modern computers appears to have reached a plateau, promoting the exploration of physics-based unconventional computational architectures beyond traditional silicon-based technology. This shift is particularly important for addressing non-deterministic polynomial-time hard (NP-hard) and NP-complete problems, which represent classes of highly complex challenges, especially in combinatorial optimization, where solution times increase exponentially with problem size, posing significant difficulties for classical computing approaches. These optimization problems [1], [2] play a crucial role in a wide range of fields, including finance [3], circuit design [4], drug discovery [5], operations [6], and scheduling [7]. Combinatorial optimization problems are known to be mappable onto ground-state search problems of the Ising model using polynomial resources [1].

The Ising Hamiltonian is given by, \[H = - \sum_{i<j}^{} J_{ij} s_i s_j - \sum_{i} h_i s_i \label{eqn:hamiltonian}\tag{1}\] where \(s_i = \pm 1\) corresponds to the \(i^{th}\) Ising spin, \(J_{ij}\) is the coupling term between spins \(s_i\) and \(s_j\), and \(h_i\) is a local bias field. The Ising machines are engineered to find configurations minimizing this Hamiltonian.

To address these computational demands of NP-hard problems, novel approaches have emerged that utilize the physical behavior of systems to perform efficient computation, such as analog Ising machines. A wide range of Ising machines have been developed using diverse physical platforms, including superconducting quantum bits [8], single-electron devices [9], nanomechanical systems [10], stochastic nanomagnets [11], CMOS circuits [12], spin waves [13]–[15], superparamagnetic tunnel junctions [16], delay-line oscillators [17], and photonic Ising machines [18], [19]. Several commercial products are also available on the market for solving combinatorial problems, including the D-Wave system, which uses superconducting technology, NTT Research’s optical Coherent Ising machine, Fujitsu’s Digital Annealer, and Hitachi’s CMOS Annealing Machine. However, each of these devices has limitations, such as the need for cryogenic temperatures, high power consumption, and confinement to laboratory environments.

Analog Ising machines using time-multiplexing have seen significant advancements, beginning with the demonstration of the optical Coherent Ising Machine (CIM) [19]–[21], this architecture overcame connectivity limitations by enabling all-to-all coupling between spins. CIMs have been successfully scaled from as few as 4 spins to over 100,000 spins. However, despite these achievements, they remain confined to laboratory settings due to their large physical footprint, high power consumption, poor temperature stability, and high cost. Some of these challenges have been partially addressed by the Spinwave Ising Machine (SWIM) [13] and surface acoustic wave-based Ising machines (SAWIM) [22], where spins are represented by spinwave or acoustic wave pulses propagating at much slower speeds, on the order of kilometers per second. However, these approaches have so far only demonstrated systems with a limited number of spins—8 in SWIM and 50 in SAWIM—and are restricted to 1-bit resolution in their coupling terms.

In this work, we introduce a time-multiplexed Ising machine based on a bulk acoustic wave delay line capable of hosting 2,048 spins. Bulk-acoustic-wave-based Ising machines (BAWIMs) offer a promising approach to building thermally stable, time-multiplexed computing systems by leveraging propagating wave packets in solid-state delay lines. Compared to state-of-the-art CIMs, BAWIM is more stable and reliable, while eliminating the need for complex frequency-stabilization systems because their thermal stability is approximately four orders of magnitude higher. We evaluate our Ising machine using arbitrary graphs from the BiqMac library and demonstrate that it achieves competitive energy results compared to state-of-the-art algorithmic approaches while effectively relaxing to low-energy states. Additionally, we report results on more complex, real-world problems, such as number partitioning and Sudoku puzzles. We show that BAWIM successfully solves the Sudoku puzzle and produces good approximate solutions for the number partitioning problem. Furthermore, we compare its performance with the simulated bifurcation algorithm and demonstrate that BAWIM outperforms it.

Figure 1: The Bulk acoustic wave based Ising machine (BAWIM). BAWIM architecture with two daisy‑chained quartz BAW delay lines for time‑multiplexed RF pulse circulation. The PSA is a parametric phase-sensitive amplifier, and the LA blocks are microwave linear amplifiers. The PSA enforces binary phases, while a measurement‑and‑feedback path reads pulse phases and injects couplings computed by the FPGA. An input and output thin-film transducer (TFT) is used to excite and receive propagating bulk acoustic waves. The inset shows the picture of the delay line. — Figure 1: **The Bulk acoustic wave based Ising machine (BAWIM).** BAWIM architecture with two daisy‑chained quartz BAW delay lines for time‑multiplexed RF pulse circulation. The PSA is a parametric phase-sensitive amplifier, and the LA blocks are microwave linear amplifiers. The PSA enforces binary phases, while a measurement‑and‑feedback path reads pulse phases and injects couplings computed by the FPGA. An input and output thin-film transducer (TFT) is used to excite and receive propagating bulk acoustic waves. The inset shows the picture of the delay line.

2 Results and discussion↩︎

2.1 Bulk acoustic waves↩︎

Bulk acoustic waves (BAWs) represent mechanical oscillations that propagate through the entirety of a material’s volume, positioning them as highly effective carriers for signal transduction in applications such as RF filtering, sensing, and delay lines within electronic systems [23], [24]. In contrast to surface acoustic waves (SAWs), which are restricted to the substrate surface and thereby exhibit heightened vulnerability to environmental contaminants, temperature variations, and fabrication defects that can impair signal integrity [25], [26], BAWs confer multiple advantages, including enhanced power handling due to volumetric energy distribution, enabling high-frequency operation often exceeding 10 GHz with reduced insertion losses and elevated Q-factors for superior frequency selectivity [24]. Furthermore, BAW devices demonstrate improved robustness against surface anomalies, better thermal stability through mechanisms like temperature-compensated designs, and greater miniaturization potential via thin-film integration [27], facilitating reliable performance in compact, high-demand scenarios such as 5G mobile communications [28] where maintaining signal fidelity in dense spectral environments is essential.

Figure 2: Experimental demonstration of BAWIM (a) Graph of the solved MAX-CUT instance with 10\% edge density, comprising 209,610 edges, each with a weight between 0 and 2^{15}; only 2,096 edges (1%) are shown for clarity. (b) Time traces of the first 200 spins during the solution run. After the 40th circulation, when the coupling is turned on, spin flips are highlighted in green. (c) Time evolution of the Ising energy across all runs, with the highlighted darker blue curve showing the average solution. During the first 40 circulations, the MFB is turned off to allow the spins to settle into random states. At the 40th circulation, the MFB is switched on, initiating the solution process. The dashed line indicates the target Ising-energy, set to 90\% of the best Ising energy obtained from simulated bifurcation. (d) Solution probabilities for the MAX-CUT problem, showing the Ising energy before the solution process begins and the final Ising energies after completion. The dashed line shows the target Ising energy that is 90\% of the simulated bifurcation’s best result. — Figure 2: **Experimental demonstration of BAWIM** **(a)** Graph of the solved MAX-CUT instance with 10\(\%\) edge density, comprising 209,610 edges, each with a weight between 0 and \(2^{15}\); only 2,096 edges (1%) are shown for clarity. **(b)** Time traces of the first 200 spins during the solution run. After the 40th circulation, when the coupling is turned on, spin flips are highlighted in green. **(c)** Time evolution of the Ising energy across all runs, with the highlighted darker blue curve showing the average solution. During the first 40 circulations, the MFB is turned off to allow the spins to settle into random states. At the 40th circulation, the MFB is switched on, initiating the solution process. The dashed line indicates the target Ising-energy, set to 90\(\%\) of the best Ising energy obtained from simulated bifurcation. **(d)** Solution probabilities for the MAX-CUT problem, showing the Ising energy before the solution process begins and the final Ising energies after completion. The dashed line shows the target Ising energy that is \(90\%\) of the simulated bifurcation’s best result.

2.2 Design of a Bulk Acoustic Wave Ising Machine↩︎

Figure 1 shows the schematic of the BAWIM, inspired by the time-multiplexing method in CIMs [19], [21], SWIMs [13], and SAWIMs [22]. At the core of the design is a quartz-based bulk acoustic wave 707 \(\mu s\) delay line, with a 3-dB bandwidth of 8 MHz centered around 20.5 MHz and with a negligible dispersion delay \(\tau _{p, disp}\) of 50 ns. Due to the resonance mechanism underlying phase-sensitive amplification, a minimum of 5–10 cycles per pulse is required, setting a lower bound on the pulse duration to approximately 250–500 ns. This constraint dominates and effectively determines the pulse period. We chose a pulse period of 665 ns with a 50\(\%\) duty cycle, corresponding to 6–7 cycles of the 20.5 MHz carrier signal per pulse. This is enforced in our design through a Mini-Circuits ZASWA-2-50DRA+ switch, shown in Fig 1. The spin states +1 and -1 are encoded in the relative phase of the pulses, with bistable phases achieved by phase-binarizing the propagating RF pulses using a phase-sensitive amplifier (PSA). The PSA comprises a phase-sensitive attenuator combined with a linear amplifier, implemented using a Mini-Circuits ZYSWA-2-50DR+ device, which behaves as a 0/1 time-gate, so the signal-component at the output is largest when the gate windows line up with the signal peaks (\(\phi = 0 ^\circ, 180 ^\circ\)) and is strongly suppressed when they line up near the zero crossings (\(\phi = \pm 90 ^\circ\)). In addition to the components mentioned above, the BAWIM uses an Analog Devices AD835 multiplier as the mixer, an Analog Devices AD8302 as the phase detector, and Mini-Circuits Z99SC-62-S+ devices as the couplers. As a single BAW delay line supports up to 1062 spins with a pulse width of 665 ns, we daisy-chain two delay lines, resulting in a total of 1.41 ms of delay, accommodating a total of 2124 spins. Of these, 2048 pulses are utilized as Ising spins, while the remaining 76 pulses are left unused.

The BAW delay line, together with linear amplifiers and the PSA, forms a multi-physics ring oscillator. The key step of the design is to establish a stable circulation of RF pulses within this loop. This requires satisfying the Barkhausen stability criteria: the loop gain must be larger than one, and the total phase shift around the loop must be an integer multiple of \(2\pi\). The delay line introduces 36 dB of attenuation, which can be compensated by the linear amplifiers, and the phase accumulation can be ensured with a proper choice of reference signal, thus satisfying the Barkhausen criteria. The PSA further narrows the stability criteria to only those signals with either phase 0 or \(\pi\) relative to the pumping signal that is twice the reference frequency.

To realize the Ising machine, we developed a measurement and feedback block (MFB) similar to the CIMs [19], [21]. Here, the portions of the circulating spin pulses are split with a 1:10 coupler, and their phase relative to a reference signal is measured by an AD8302 phase detector. The resulting signal is digitized by an 8-bit AD9280 ADC and processed by an AMD Zynq ZCU104 FPGA. The FPGA output is then converted back to analog using an AD9708 DAC, multiplied with the reference signal using an AD835 chip, and reinjected into the corresponding spin via a coupler.

The FPGA calculates the coupling pulses (\(c_i\)) based on the coupling matrix \(J_{ij}\), and spins \(s_j\),

\[c_i = \sum_{j=1}^{N} J_{ij} s_j\] Here the coupling terms \(J_{ij}\) can have 15-bit resolution (0 to \(+2^{15}\)). The \(c_i\) controls the phase and amplitude of the signal being injected to the circulating RF signal. The amplitude of the coupling signal can be tuned using the potentiometer of the DAC, allowing an additional control to ensure the injected signal remains within \(5-30\%\) of the circulating RF signal and does not overpower it, avoiding chaotization of the oscillatory circuit.

2.3 Performance of the Ising Machine↩︎

We evaluate the BAWIM performance using MAX-CUT problems on graph instances from the BiqMac library [29], [30].Figure 2 (a) shows a 2048-node graph with an edge density of 10\(\%\) having 209,610 edges, where each edge can have a weight ranging from 0 to \(2^{15}\). The time evolution of the first 200 spins is presented in Fig. 2 (b), where it is observed that the BAWIM immediately begins searching for lower energy states and gradually slows down as it reaches better solutions. Figure 2 (c) shows the temporal evolution of the Ising energy for multiple runs; the similar energy evolution over different runs demonstrates the reliability of our design, and it also indicates the target Ising energy based on the simulated bifurcation algorithm (details in the next section). For consecutive runs, we turn off and on the amplification in the loop with a \(3\%\) duty cycle control signal with a repetition period of 2 seconds, which allows the BAW pulse to decay completely. The quantitative results of the MAX-CUT problem are illustrated in Fig. 2 (d) by a histogram of the Ising energies and MAX-CUT scores for over 100 runs, along with the target Ising energy. The performance of BAWIM depends on the coupling strength, and there exists an optimal coupling strength depending on the problem instance. For a MAX-CUT problem, the number of cuts is given by \[Cuts = -\frac{1}{2} \sum_{i<j}^{} J_{ij} - \frac{1}{2} H,\] where H is the Hamiltonian (Ising energy) in Eq. 1 .

Figure 3: Comparison with heated-ballistic simulated bifurcation algorithm (a) Time evolution of the Ising energy for BAWIM, HbSB, and the projected performance of BAWIM with a high central frequency delay line, with a target Ising energy of 90 \% of HbSB’s best result. (b) Success rate and time-to-target (TTT) of BAWIM relative to the best MAX-CUT score of HbSB. The blue curve shows the success rate computed from n = 110 independent BAWIM runs. The orange curve shows the mean TTT computed over successful runs, and the error bars indicate \pm1 standard deviation. — Figure 3: **Comparison with heated-ballistic simulated bifurcation algorithm** **(a)** Time evolution of the Ising energy for BAWIM, HbSB, and the projected performance of BAWIM with a high central frequency delay line, with a target Ising energy of \(90 \%\) of HbSB’s best result. **(b)** Success rate and time-to-target (TTT) of BAWIM relative to the best MAX-CUT score of HbSB. The blue curve shows the success rate computed from n = 110 independent BAWIM runs. The orange curve shows the mean TTT computed over successful runs, and the error bars indicate \(\pm1\) standard deviation.

Comparison with Simulated bifurcation algorithm↩︎

In order to perform comparative benchmarking of the BAWIM, we use the Simulated Bifurcation (SB) algorithm [31], [32], which has outperformed CIMs and simulated annealing in both solution quality and time-to-solution [31]. The SB algorithm has been further enhanced in subsequent works [33], [34] by incorporating principles from classical mechanics and thermal heating. In this study, we use the Heated Ballistic Simulated Bifurcation (HbSB) variant of SB, as it has demonstrated improved performance among various versions of SB [34]. Using HbSB, we conducted 1,000 simulations to determine the minimum Ising energy and the maximum MAX-CUT score, as summarized in Table 1. Figure 3 (a) compares the time evolution of the best solution obtained by HbSB and BAWIM. Our machine required 341 ms to reach 90\(\%\) of HbSB’s best Ising energy, whereas HbSB achieved the same energy in 102 ms. The HbSB simulations were run on a modern system with an AMD Ryzen 7 9700X 8-core (5.30 GHz) CPU and 64 GB of RAM. The performance of BAWIM can be significantly enhanced by utilizing a BAW delay line with a higher central frequency. In our current setup, we use a delay line with a 20.5 MHz central frequency, whereas BAW delay lines can operate at much higher frequencies up to 26 GHz. A catalog [35] from Teledyne lists the publicly available devices, among them, the MBJ-1018 device, operating at 16.455 GHz, can accommodate a similar number of spins to our present device. Operating the design at this higher frequency may require custom RF components as well as an FPGA capable of operating at higher frequencies. Figure 3 (a) shows the projected runtime improvement using this 16.455 GHz delay line, where the target energy is estimated to be reached in just 0.462 ms, which is about 750 times faster than the current implementation.

Figure 3 (b) shows the success rate and time-to-target, relative to the best MAX-CUT score obtained by HbSB. Across all runs, BAWIM consistently achieves 99\(\%\) of HbSB’s optimal MAX-CUT score, with longer times required as the target increases. The BAWIM performance also varies with the problem size.

Dependence on the edge density↩︎

To comprehensively benchmark our system, we evaluate its performance on problems with varying edge densities. Specifically, we test BAWIM using problems with edge densities of \(25\%\), \(50\%\), \(75\%\), and \(100\%\), generated from the BiqMac library. Figure 4 presents the quantitative results based on over 100 solutions for each problem, showcasing the MAX-CUT scores in comparison to the HbSB’s best result. We also evaluated these problems using the HbSB algorithm. Table 1 reports the best observed Ising energy and MAX-CUT score from 1,000 simulations of the HbSB algorithm; it also includes the best performance achieved by our system among 100 runs. Here, we observe that BAWIM attains up to \(95\%\) of HbSB’s best performance in terms of Ising energy, and up to \(99.99\%\) in terms of the MAX-CUT score.

Figure 4 (a) illustrates the time-to-target for problems with varying edge densities, using a target set at \(90\%\) of HbSB’s best result. For the \(10\%\) edge density problem, the time-to-target is approximately 341 ms and reaches the target \(68\%\) of the time, whereas the fully connected (\(100\%\) edge density) problem requires around 774 ms and achieves a success rate of \(97\%\).

Figure 4: (a) Time to solution and success rate for different edge density problems, with target Ising energy of 90\% of HbSB’s best result. The orange curve shows the success rate computed from n = 110 independent BAWIM runs. The blue curve shows the mean TTT computed over successful runs, and the error bars indicate \pm1 standard deviation. Solution probabilities for MAX-CUT problems with various edge densities (b) 25\% edge density. (c) 50\% edge density. (d) 75\% edge density. (e) 100\% edge density. — Figure 4: **(a) Time to solution and success rate for different edge density problems, with target Ising energy of \(90\%\) of HbSB’s best result.** The orange curve shows the success rate computed from n = 110 independent BAWIM runs. The blue curve shows the mean TTT computed over successful runs, and the error bars indicate \(\pm1\) standard deviation. **Solution probabilities for MAX-CUT problems with various edge densities** **(b)** 25\(\%\) edge density. **(c)** 50\(\%\) edge density. **(d)** 75\(\%\) edge density. **(e)** 100\(\%\) edge density.

Table 1: Comparison of BAWIM and HbSB algorithms in terms of Ising energy across problems with different edge densities. The BAWIM Ising energies are expressed as percentages of the corresponding HbSB values.
Edge Density	Best Ising Energy
Edge Density	HbSB	BAWIM
10\(\%\)	-803,331,776	\(93.10 \%\)
25\(\%\)	-1,201,604,480	\(93.59 \%\)
50\(\%\)	-1,491,890,816	\(95.52 \%\)
75\(\%\)	-1,535,851,136	\(95.74 \%\)
100\(\%\)	-1,339,657,856	\(94.52 \%\)

2.4 Number Partitioning↩︎

In this section, we focus on the number partitioning problem (NPP), which is one of Karp’s 21 NP-complete problems [36]. It is defined as the task of dividing a given set of positive integers into two subsets such that the absolute difference between their sums is minimized, an example is shown in Fig. 5 (a). The NPP has practical applications ranging from multiprocessor scheduling and minimizing the size and delay of the VLSI circuits [37], [38], to public-key cryptography [39], and in choosing sides in a ball game [40]. Furthermore, both the exact cover and the knapsack problems can be reduced linearly [36], [41] to NPP without any constraints. This indicates that solving the NPP would also address these other NP-hard problems.

The NPP can be described as, \[E(\mathcal{A},\mathcal{B}) = \left| \sum_{i \in \mathcal{A}} a_i - \sum_{i \in \mathcal{B}} a_i \right|\] where \(a_1, a_2,....a_N\) are the list of positive integers, and \(\mathcal{A}\), \(\mathcal{B}\) are the subsets of this list, and E represents the difference between the two groups, with E=0(1) representing the perfect partition when \(\sum a_i\) is even(odd).

Let \(s=+1\) represent subset \(\mathcal{A}\) and \(s=-1\) represent subset \(\mathcal{B}\), then the above equation be written as \[E(s) = \left| \sum_{i=1}^{N} a_i s_i \right|\]
So, the coupling terms for the Ising model for the NPP can be written as \(J_{ij}=a_i a_j\) with bias field \(h=0\). This Hamiltonian for NPP is similar to the Mattis spin glass Hamiltonian [42], [43].

To test the NPP, we generated seven distinct integer sets with sizes ranging from 32 to 2048, where each integer lies within the range of 1 to 178 and was generated randomly in Python. In BAWIM, the coupling elements \(J_{ij}\) are limited to a 15-bit resolution, which limits the maximum integer values that can be used in the NPP. These problem instances were evaluated using both the HbSB algorithm and the BAWIM. The results are summarized in Fig. 5 (b), here each set was tested for 10,000 runs with the HbSB algorithm and around 130 to 700 runs with the BAWIM. In BAWIM, smaller-node problems were run in parallel as replicas to get more quantitative results. In BAWIM, each solution was executed for 2 seconds, and the smallest difference obtained during each run was reported. Here, the exact solution corresponds to a difference of 0, while the \(99.9\%\) approximate solution indicates that the difference is less than \(0.1\%\) of the maximum possible value (i.e., the total sum of all numbers). In all cases, BAWIM outperforms the HbSB algorithm, achieving a \(100\%\) success rate in most instances when considering the \(99.9\%\) approximate solution.

Figure 5: Number partitioning problem (a) Example of a number partitioning problem, represented as an all-to-all connected graph, where each spin denotes the assignment of a number to one of two subsets with equal sums. (b) The success rates of BAWIM and the HbSB algorithm are shown across problem sizes ranging from 32 to 2048, evaluated for different levels of approximate solution. — Figure 5: **Number partitioning problem** **(a)** Example of a number partitioning problem, represented as an all-to-all connected graph, where each spin denotes the assignment of a number to one of two subsets with equal sums. **(b)** The success rates of BAWIM and the HbSB algorithm are shown across problem sizes ranging from 32 to 2048, evaluated for different levels of approximate solution.

2.5 Sudoku↩︎

Sudoku is a famous puzzle game played on a \(9 \times 9\) grid, where each cell can contain a number from 1 to 9. At the start, some cells are pre-filled with “clues,” while the rest remain blank. The goal is to complete the grid so that (i) every row, (ii) every column, and (iii) every 3 × 3 sub-grid (“block”) each contains all numbers from 1 to 9 exactly once. The Sudoku puzzle originated in the 1970s as “Latin squares”, later adapted into its modern form with the name “Number Place”, and eventually rose to global popularity as “Sudoku” [44], [45]. Beyond attracting puzzle enthusiasts, it has also become a frequent subject of study for mathematicians and computer scientists [44]–[49], likely owing to its simplicity and accessibility to a non-technical audience. For generalized grid sizes of \(n^2 \times n^2\), the problem has been proven to be NP-complete [50].

To adapt the Sudoku problem to the Ising model, we represent the digits 1 through 9 in each cell using a one-hot encoding scheme, since the Ising spins are inherently binary. This requires 9 spins per cell, leading to a total of \(81 \times 9 = 729\) spins for the entire grid. This encoding adds a new rule that only one number should be present in each cell. The coupling matrix is constructed by extending the one-hot encoding to incorporate the Sudoku rules, as outlined in [51]. Clues are embedded in the bias field \(h\), for a cell with a given clue, the bias of the corresponding spin is lowered while the biases of the other eight spins are raised, resulting in the clue being encouraged.

Figure 6 (a) illustrates the evolution of the Ising energy and the number of rule-violating cells during the Sudoku-solving process using BAWIM. Figures 6 (b) and (c) display intermediate puzzle states corresponding to 96.71 \(\%\) and 99.42 \(\%\) of the ground-state energy, respectively. Notably, even when the Ising energy approaches the ground state, the Sudoku solution remains far from correct. The fully solved puzzle is shown in Fig. 6 (d), where black numbers represent the given clues and green numbers indicate the completed solution. We also tested the Sudoku problems using the HbSB algorithm; however, it failed to produce correct solutions.The Sudoku benchmark is included primarily to demonstrate the programmability of BAWIM for dense, non-native constraint-satisfaction problems. This highlights that, unlike other optimization problems where near-optimal states can still be useful, Sudoku requires exact constraint satisfaction, therefore, a near-ground-state Ising energy may still correspond to an invalid solution.

Figure 6: Sudoku puzzle (a) Evolution of the Ising energy and the number of cells with rule violations over time, the inset shows a zoomed-in view. The three vertical lines indicate the moments when the Ising energy reaches 96.71\%, 99.42\%, and 100\% of the ground state. The corresponding Sudoku configurations are shown in (b), (c), and (d). The background colors highlight different types of issues associated with the cell values. The black numbers denote the original clues, multiple numbers within a single cell indicate multiple candidates violation, and the green numbers show the filled-in solutions. — Figure 6: **Sudoku puzzle** **(a)** Evolution of the Ising energy and the number of cells with rule violations over time, the inset shows a zoomed-in view. The three vertical lines indicate the moments when the Ising energy reaches 96.71\(\%\), 99.42\(\%\), and 100\(\%\) of the ground state. The corresponding Sudoku configurations are shown in **(b)**, **(c)**, and **(d)**. The background colors highlight different types of issues associated with the cell values. The black numbers denote the original clues, multiple numbers within a single cell indicate multiple candidates violation, and the green numbers show the filled-in solutions.

2.6 Estimation of Power consumption↩︎

The total power consumption of BAWIM is approximately 9.57 W, primarily dominated by the FPGA and the five ZFL-500LN+ amplifiers. The AMD Zynq ZCU104 FPGA consumes about 5.2 W, with 2.40 W attributed to the programmable logic and 2.80 W to the processing system. The five ZFL-500LN+ amplifiers together draw 3.75 W, replacing them with specialized low-power amplifiers could reduce this to 50–100 mW. The remaining components, including the phase detector, RF switches, ADC, DAC, and multiplier, collectively contribute about 0.6 W. By utilizing low-power amplifiers, the power consumption of BAWIM can be reduced to 5.9 W.

2.7 Comparison with CIMs↩︎

The BAWIM builds on CIM foundations with major advances in functionality and architecture. Prior CIM implementations faced some limitations that our system avoids, for example, a 2000-spin CIM [21] required extensive phase-locking and stabilization, restricting usable spins to 2048 out of 5082 DOPOs, with the rest devoted to auxiliary roles. In contrast, BAWIM uses 2048 of 2124 RF pulses directly as spins, with the remaining 76 left unused. Similarly, 100,000-spin CIM [19] also suffered from fiber-cavity phase instability, inconsistent solution times, and heavy post-selection, discarding up to half of the runs due to incomplete phase erasure. Our system operates stably without post-selection, improving both efficiency and throughput. Moreover, while CIM performance depended strongly on pump scheduling, BAWIM does not use such scheduling, yielding more predictable behavior.

Due to the disparity in carrier frequencies, BAWIM exhibits significantly higher efficiency in terms of carrier periods per spin and spin duty cycle. For comparison, the first CIM involved approximately \(1\times 10^6\) carrier periods, while the more recent 100,000-spin CIM used around 6,000 with a \(15\%\) duty cycle. In contrast, BAWIM operates with just 6 carrier periods per spin and with a \(50\%\) duty cycle.

In addition to the differences mentioned above, a significant challenge with CIMs is their sensitivity to temperature. The 100,000-spin CIM [19] has a temperature coefficient of phase accumulation of \(1.73 \times 10^7\) \(deg/{^\circ} C\), meaning that a temperature change of just \(1^\circ C\) can cause a shift in group delay of approximately 250 ps in a 5-km fiber delay line. To address this, CIMs must operate within a tightly controlled thermal environment with an accuracy of \(\pm 0.05^\circ C\). The authors achieve this by placing the fiber inside a container with thick, hollow walls filled with water. Multiple thermistors and Peltier devices are used for active temperature regulation, and the setup is insulated with styrofoam and enclosed in an aluminum box. In contrast, our BAWIM exhibits a temperature coefficient of phase accumulation of only 780 \({deg /^\circ}C\), which is \(2.21 \times 10^4\) times lower than that of the CIM.As a result, it requires no environmental or temperature control and can operate reliably at room temperature on a standard tabletop setup. Its ease of use, thermal stability, low power consumption, benchtop design, and the use of only off-the-shelf RF components make BAWIM a more advantageous and commercially viable option for combinatorial solvers.

3 Conclusion↩︎

We presented the implementation of a time-multiplexed Ising machine based on bulk acoustic waves propagating through a quartz-based solid-state delay line. The BAWIM system features 2048 spins with 15-bit coupling resolution and with a circulation time of 1.41 ms. We demonstrated its performance by solving the MAX-CUT problem on arbitrary graphs, achieving a solution time of 341 ms with the potential to reduce this to 0.46 ms. The BAWIM was benchmarked against graphs of varying densities from the BiqMac library. Beyond MAX-CUT, we solve more complex tasks such as number partitioning and Sudoku, highlighting its practical utility for real-world problems. Our design offers orders-of-magnitude higher thermal stability compared to state-of-the-art coherent Ising machines and surpasses the simulated bifurcation algorithm’s performance in tackling complex problems. Overall, the BAWIM provides a reliable, benchtop, thermally stable, and energy-efficient architecture with strong potential for further improvements in solution time and power consumption, with a promising prospect for commercialization.

4 Methods↩︎

Sample: The BAW delay medium is a fused-quartz disc shaped as a tetradecagon, measuring 9.5 cm across opposite sides. The microwave-to-acoustic wave transducers are made of Y-cut quartz crystal with silver electrodes.

Electrical measurements: The frequency and group delay of the quartz delay line are derived from S-parameters measured using a vector network analyzer. Measurements were performed over a frequency range of 10–30 MHz with a resolution of 0.2 kHz. Group delay is determined by taking the negative derivative of the phase shift with respect to frequency, and both the loss and time-delay curves are smoothed using a moving average filter with a window of 4000 points. The phase-sensitive amplifier was characterized by obtaining the envelope of an amplified signal with constant amplitude and a frequency detuned by 1 kHz from the reference signal. This detuning introduced a continuous phase drift with a 1 ms period, enabling the extraction of a smooth phase-sensitive amplification curve.

Acknowledgments↩︎

This work was supported by a Knut and Alice Wallenberg Foundation WALP grant, KAW 253129326, a Horizon 2020 research and innovation program ERC advanced grant no. 835068 “TOPSPIN”, an ERC proof of concept grant no. 101069424 “SPINTOP”, and the Marie Skłodowska-Curie grant agreement no. 101111429 “SWIM”.

Author Contributions↩︎

A.L. and J.Å. conceived the concept. V.V. and A.L. designed the circuit. V.V performed the measurements and analyzed the data. V.V., R.O., V.G., and R.K. performed theoretical calculations. J.Å. managed the project. All co-authors contributed to the manuscript, the discussion, and the analysis of the results.

References↩︎

[1]

F. Barahona, “On the computational complexity of ising spin glass models,” Journal of Physics A: Mathematical and General, vol. 15, no. 10, p. 3241, 1982.

[2]

D.-Z. Du and P. M. Pardalos, Eds., Handbook of combinatorial optimization, vol. 4. New York, NY: Springer, 1998.

[3]

O. H. Ibarra and C. E. Kim, “Fast approximation algorithms for the knapsack and sum of subset problems,” Journal of the ACM (JACM), vol. 22, no. 4, pp. 463–468, 1975.

[4]

F. Barahona, M. Grötschel, M. Jünger, and G. Reinelt, “An application of combinatorial optimization to statistical physics and circuit layout design,” Operations Research, vol. 36, no. 3, pp. 493–513, 1988.

[5]

D. J. Earl and M. W. Deem, “Parallel tempering: Theory, applications, and new perspectives,” Physical Chemistry Chemical Physics, vol. 7, no. 23, pp. 3910–3916, 2005.

[6]

V. Černỳ, “Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm,” Journal of optimization theory and applications, vol. 45, no. 1, pp. 41–51, 1985.

[7]

E. K. Burke, P. De Causmaecker, G. V. Berghe, and H. Van Landeghem, “The state of the art of nurse rostering,” Journal of scheduling, vol. 7, no. 6, pp. 441–499, 2004.

[8]

M. W. Johnson et al., “Quantum annealing with manufactured spins,” Nature, vol. 473, no. 7346, pp. 194–198, 2011.

[9]

K. Nishiguchi and A. Fujiwara, “Single-electron circuit for stochastic data processing using nano-MOSFETs,” in 2007 IEEE international electron devices meeting, 2007, pp. 791–794.

[10]

I. Mahboob, H. Okamoto, and H. Yamaguchi, “An electromechanical ising hamiltonian,” Science advances, vol. 2, no. 6, p. e1600236, 2016.

[11]

B. Sutton, K. Y. Camsari, B. Behin-Aein, and S. Datta, “Intrinsic optimization using stochastic nanomagnets,” Scientific reports, vol. 7, no. 1, p. 44370, 2017.

[12]

W. Whitehead, Z. Nelson, K. Y. Camsari, and L. Theogarajan, “CMOS-compatible ising and potts annealing using single-photon avalanche diodes,” Nature Electronics, vol. 6, no. 12, pp. 1009–1019, 2023.

[13]

A. Litvinenko et al., “A spinwave ising machine,” Communications Physics, vol. 6, no. 1, p. 227, 2023.

[14]

V. H. González, A. Litvinenko, R. Khymyn, and J. Åkerman, “Global biasing using a hardware-based artificial zeeman term in spinwave ising machines,” Applied Physics Letters, vol. 124, no. 9, 2024.

[15]

V. H. González, A. Litvinenko, A. Kumar, R. Khymyn, and J. Åkerman, “Spintronic devices as next-generation computation accelerators,” Current Opinion in Solid State and Materials Science, vol. 31, p. 101173, 2024.

[16]

J. Si et al., “Energy-efficient superparamagnetic ising machine and its application to traveling salesman problems,” Nature Communications, vol. 15, no. 1, p. 3457, 2024.

[17]

R. V. Ovcharov, V. H. González, A. Litvinenko, J. Åkerman, and R. S. Khymyn, “A numerical model for time-multiplexed ising machines based on delay-line oscillators,” arXiv preprint arXiv:2406.07197, 2024.

[18]

Y. Yamamoto et al., “Coherent ising machines—optical neural networks operating at the quantum limit,” npj Quantum Information, vol. 3, no. 1, p. 49, 2017.

[19]

T. Honjo et al., “100,000-spin coherent ising machine,” Science advances, vol. 7, no. 40, p. eabh0952, 2021.

[20]

H. Takesue et al., “Finding independent sets in large-scale graphs with a coherent ising machine,” Science Advances, vol. 11, no. 7, p. eads7223, 2025.

[21]

T. Inagaki et al., “A coherent ising machine for 2000-node optimization problems,” Science, vol. 354, no. 6312, pp. 603–606, 2016.

[22]

A. Litvinenko, R. Khymyn, R. Ovcharov, and J. Åkerman, “A 50-spin surface acoustic wave ising machine,” Communications Physics, vol. 8, no. 1, pp. 1–11, 2025.

[23]

K. Hashimoto, Ed., RF bulk acoustic wave filters for communications. Norwood, MA: Artech House, 2009.

[24]

Y. Liu, Y. Cai, Y. Zhang, A. Tovstopyat, S. Liu, and C. Sun, “Materials, design, and characteristics of bulk acoustic wave resonator: A review,” Micromachines, vol. 11, no. 7, p. 630, 2020.

[25]

D. Mandal and S. Banerjee, “Surface acoustic wave (SAW) sensors: Physics, materials, and applications,” Sensors, vol. 22, no. 3, p. 820, 2022.

[26]

Z. Tang et al., “A review of surface acoustic wave sensors: Mechanisms, stability and future prospects,” Sensor Review, vol. 44, no. 3, pp. 249–266, 2024.

[27]

M. D. Terzieva, “Overview of bulk acoustic wave technology and its applications.” Electrotechnica & Electronica (E+ E), vol. 51, 2016.

[28]

R. Aigner, G. Fattinger, M. Schaefer, K. Karnati, R. Rothemund, and F. Dumont, “BAW filters for 5G bands,” in 2018 IEEE international electron devices meeting (IEDM), 2018, pp. 14–5.

[29]

G. Rinaldi, “Rudy graph generator.” www-user.tu-chemnitz.de/˜helmberg/rudy.tar.gz, 1996.

[30]

C. Helmberg and F. Rendl, “A spectral bundle method for semidefinite programming,” SIAM Journal on Optimization, vol. 10, no. 3, pp. 673–696, 2000.

[31]

H. Goto, K. Tatsumura, and A. R. Dixon, “Combinatorial optimization by simulating adiabatic bifurcations in nonlinear hamiltonian systems,” Science advances, vol. 5, no. 4, p. eaav2372, 2019.

[32]

R. Ageron, T. Bouquet, and L. Pugliese, “Simulated bifurcation (SB) algorithm for python.” https://github.com/bqth29/simulated-bifurcation-algorithm, 2025.

[33]

H. Goto et al., “High-performance combinatorial optimization based on classical mechanics,” Science Advances, vol. 7, no. 6, p. eabe7953, 2021.

[34]

T. Kanao and H. Goto, “Simulated bifurcation assisted by thermal fluctuation,” Communications Physics, vol. 5, no. 1, p. 153, 2022.

[35]

Teledyne, “BAW devices product selection guide.” https://www.teledynedefenseelectronics.com/wireless/Wireless Brochures/BAW Devices Product Selection Guide.pdf.

[36]

R. M. Karp, “Reducibility among combinatorial problems,” in Complexity of computer computations: Proceedings of a symposium on the complexity of computer computations, held march 20–22, 1972, at the IBM thomas j. Watson research center, yorktown heights, new york, and sponsored by the office of naval research, mathematics program, IBM world trade corporation, and the IBM research mathematical sciences department, R. E. Miller, J. W. Thatcher, and J. D. Bohlinger, Eds. Boston, MA: Springer US, 1972, pp. 85–103.

[37]

Jr. Coffman Edward Grady and G. S. Lueker, Probabilistic analysis of packing and partitioning algorithms. New York; Chichester: WileyInterscience, 1991.

[38]

L.-H. Tsai, “Asymptotic analysis of an algorithm for balanced parallel processor scheduling,” SIAM Journal on Computing, vol. 21, no. 1, pp. 59–64, 1992.

[39]

R. Merkle and M. Hellman, “Hiding information and signatures in trapdoor knapsacks,” IEEE transactions on Information Theory, vol. 24, no. 5, pp. 525–530, 1978.

[40]

B. Hayes, “Computing science: The easiest hard problem,” American Scientist, vol. 90, no. 2, pp. 113–117, 2002.

[41]

Z. Li, T. Seidel, D. Leib, M. Bortz, and R. Heese, “Efficient solution of the number partitioning problem on a quantum annealer: A hybrid quantum-classical decomposition approach,” Journal of Heuristics, vol. 31, no. 2, p. 21, 2025.

[42]

D. Mattis, “Solvable spin systems with random interactions,” Physics Letters A, vol. 56, no. 5, pp. 421–422, 1976.

[43]

T. Graß, D. Raventós, B. Juliá-Dı́az, C. Gogolin, and M. Lewenstein, “Quantum annealing for the number-partitioning problem using a tunable spin glass of ions,” Nature communications, vol. 7, no. 1, p. 11524, 2016.

[44]

T. Resnick, “Sudoku at the intersection of classical and quantum computing,” Department of Computer Science, The University of Auckland, New Zealand, 2014.

[45]

S. Mücke, “A simple QUBO formulation of sudoku,” in Proceedings of the genetic and evolutionary computation conference companion, 2024, pp. 1958–1962.

[46]

R. Pignari, V. Fra, E. Macii, and G. Urgese, “Efficient solution validation of constraint satisfaction problems on neuromorphic hardware: The case of sudoku puzzles,” IEEE Transactions on Artificial Intelligence, 2025.

[47]

A. Shukla, M. Erementchouk, and P. Mazumder, “Non-binary dynamical ising machines for combinatorial optimization,” Physica D: Nonlinear Phenomena, p. 134809, 2025.

[48]

M. Ercsey-Ravasz and Z. Toroczkai, “The chaos within sudoku,” Scientific reports, vol. 2, no. 1, pp. 1–8, 2012.

[49]

I. Lynce and J. Ouaknine, “Sudoku as a SAT problem.” in AI&m, 2006.

[50]

T. Yato and T. Seta, “Complexity and completeness of finding another solution and its application to puzzles,” IEICE transactions on fundamentals of electronics, communications and computer sciences, vol. 86, no. 5, pp. 1052–1060, 2003.

[51]

Y. Minato, “Solving the world’s most difficult SUDOKU problem using ising model on javascript — minatoyuichiro.medium.com.” https://minatoyuichiro.medium.com/sudoku-solving-the-worlds-most-difficult-sudoku-problem-using-ising-model-on-javascript-c9b5add0e5c0.

A 2048-spin bulk acoustic wave Ising machine for number partitioning and Sudoku