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

# A 2048-spin bulk acoustic wave Ising machine for number partitioning and Sudoku
Source: [https://arxiv.org/html/2607.02112](https://arxiv.org/html/2607.02112)
\[1\]Venkatesh Vadde \[1,2\]Artem Litvinenko \[1,3,4\]Johan Åkerman

1\]Department of Physics, University of Gothenburg, Gothenburg, Sweden 2\]Department of Physics, Oakland University, Rochester, Michigan, USA 3\]Center for Science and Innovation in Spintronics, Tohoku University, Sendai, Japan 4\]Research Institute of Electrical Communication, Tohoku University, Sendai, Japan

###### Abstract

Optical coherent Ising machines based on time\-multiplexing have demonstrated significant progress in terms of connectivity and spin scalability\. However, they are constrained by large physical footprints, high power consumption, poor thermal stability, and high cost\. Here, we present a time\-multiplexed Ising machine leveraging propagating wave packets in solid\-state delay lines at microwave frequencies, enabling thermally stable, robust, low\-power, tabletop, and affordable design\. We use two serially connected 20\.5 MHz, 707μ​s\\mu sbulk acoustic wave delay lines supporting 2,048 spins\. Our design provides all\-to\-all connectivity with 15\-bit coupling resolution and finds approximate MAX\-CUT solutions in 341 ms, potentially scalable to sub\-ms by using higher frequency delay lines\. Additionally, we demonstrate solutions to number partitioning and Sudoku problems\. Compared with state\-of\-the\-art Coherent Ising machines, our machine exhibits four orders of magnitude higher thermal stability\. Against the simulated bifurcation algorithm, our design achieves comparable results on the MAX\-CUT problem, while outperforming it on the more complex number\-partitioning and Sudoku problems\.

## 1Introduction

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\[barahona1982computational,du1998handbook\]play a crucial role in a wide range of fields, including finance\[ibarra1975fast\], circuit design\[barahona1988application\], drug discovery\[earl2005parallel\], operations\[vcerny1985thermodynamical\], and scheduling\[burke2004state\]\. Combinatorial optimization problems are known to be mappable onto ground\-state search problems of the Ising model using polynomial resources\[barahona1982computational\]\.

The Ising Hamiltonian is given by,

H=−∑i<jJi​j​si​sj−∑ihi​siH=\-\\sum\_\{i<j\}J\_\{ij\}s\_\{i\}s\_\{j\}\-\\sum\_\{i\}h\_\{i\}s\_\{i\}\(1\)wheresi=±1s\_\{i\}=\\pm 1corresponds to theit​hi^\{th\}Ising spin,Ji​jJ\_\{ij\}is the coupling term between spinssis\_\{i\}andsjs\_\{j\}, andhih\_\{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\[johnson2011quantum\], single\-electron devices\[nishiguchi2007single\], nanomechanical systems\[mahboob2016electromechanical\], stochastic nanomagnets\[sutton2017intrinsic\], CMOS circuits\[whitehead2023cmos\], spin waves\[litvinenko2023spinwave,gonzalez2024global,gonzalez2024spintronic\], superparamagnetic tunnel junctions\[si2024energy\], delay\-line oscillators\[ovcharov2024numerical\], and photonic Ising machines\[yamamoto2017coherent,honjo2021100\]\. 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\)\[takesue2025finding,honjo2021100,inagaki2016coherent\], 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\)\[litvinenko2023spinwave\]and surface acoustic wave\-based Ising machines \(SAWIM\)\[litvinenko202550\], 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\.

![Refer to caption](https://arxiv.org/html/2607.02112v1/Figure1.png)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\.
## 2Results and discussion

### 2\.1Bulk 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\[hashimoto2009rf,liu2020materials\]\. 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\[mandal2022surface,tang2024review\], 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\[liu2020materials\]\. 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\[terzieva2016overview\], facilitating reliable performance in compact, high\-demand scenarios such as 5G mobile communications\[aigner2018baw\]where maintaining signal fidelity in dense spectral environments is essential\.

![Refer to caption](https://arxiv.org/html/2607.02112v1/Figure2.png)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 and2152^\{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 is90%90\\%of the simulated bifurcation’s best result\.
### 2\.2Design of a Bulk Acoustic Wave Ising Machine

Figure[1](https://arxiv.org/html/2607.02112#S1.F1)shows the schematic of the BAWIM, inspired by the time\-multiplexing method in CIMs\[honjo2021100,inagaki2016coherent\], SWIMs\[litvinenko2023spinwave\], and SAWIMs\[litvinenko202550\]\. At the core of the design is a quartz\-based bulk acoustic wave 707μ​s\\mu sdelay line, with a 3\-dB bandwidth of 8 MHz centered around 20\.5 MHz and with a negligible dispersion delayτp,d​i​s​p\\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](https://arxiv.org/html/2607.02112#S1.F1)\. 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 \(ϕ=0∘,180∘\\phi=0^\{\\circ\},180^\{\\circ\}\) and is strongly suppressed when they line up near the zero crossings \(ϕ=±90∘\\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 of2​π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π\\pirelative 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\[honjo2021100,inagaki2016coherent\]\. 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 \(cic\_\{i\}\) based on the coupling matrixJi​jJ\_\{ij\}, and spinssjs\_\{j\},

ci=∑j=1NJi​j​sjc\_\{i\}=\\sum\_\{j=1\}^\{N\}J\_\{ij\}s\_\{j\}\(2\)Here the coupling termsJi​jJ\_\{ij\}can have 15\-bit resolution \(0 to\+215\+2^\{15\}\)\. Thecic\_\{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 within5−30%5\-30\\%of the circulating RF signal and does not overpower it, avoiding chaotization of the oscillatory circuit\.

### 2\.3Performance of the Ising Machine

We evaluate the BAWIM performance using MAX\-CUT problems on graph instances from the BiqMac library\[rudy\_gen,helmberg2000spectral\]\.Figure[2](https://arxiv.org/html/2607.02112#S2.F2)\(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 to2152^\{15\}\. The time evolution of the first 200 spins is presented in Fig\.[2](https://arxiv.org/html/2607.02112#S2.F2)\(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](https://arxiv.org/html/2607.02112#S2.F2)\(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 a3%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](https://arxiv.org/html/2607.02112#S2.F2)\(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

C​u​t​s=−12​∑i<jJi​j−12​H,Cuts=\-\\frac\{1\}\{2\}\\sum\_\{i<j\}J\_\{ij\}\-\\frac\{1\}\{2\}H,\(3\)where H is the Hamiltonian \(Ising energy\) in Eq\.[1](https://arxiv.org/html/2607.02112#S1.E1)\.

![Refer to caption](https://arxiv.org/html/2607.02112v1/Figure3.png)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 of90%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±1\\pm 1standard deviation\.#### Comparison with Simulated bifurcation algorithm

In order to perform comparative benchmarking of the BAWIM, we use the Simulated Bifurcation \(SB\) algorithm\[goto2019combinatorial,Ageron\_Simulated\_Bifurcation\_SB\_2023\], which has outperformed CIMs and simulated annealing in both solution quality and time\-to\-solution\[goto2019combinatorial\]\. The SB algorithm has been further enhanced in subsequent works\[goto2021high,kanao2022simulated\]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\[kanao2022simulated\]\. Using HbSB, we conducted 1,000 simulations to determine the minimum Ising energy and the maximum MAX\-CUT score, as summarized in Table[1](https://arxiv.org/html/2607.02112#S2.T1)\. Figure[3](https://arxiv.org/html/2607.02112#S2.F3)\(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\[teledyne\_BAW\]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](https://arxiv.org/html/2607.02112#S2.F3)\(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](https://arxiv.org/html/2607.02112#S2.F3)\(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 of25%25\\%,50%50\\%,75%75\\%, and100%100\\%, generated from the BiqMac library\. Figure[4](https://arxiv.org/html/2607.02112#S2.F4)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](https://arxiv.org/html/2607.02112#S2.T1)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 to95%95\\%of HbSB’s best performance in terms of Ising energy, and up to99\.99%99\.99\\%in terms of the MAX\-CUT score\.

Figure[4](https://arxiv.org/html/2607.02112#S2.F4)\(a\) illustrates the time\-to\-target for problems with varying edge densities, using a target set at90%90\\%of HbSB’s best result\. For the10%10\\%edge density problem, the time\-to\-target is approximately 341 ms and reaches the target68%68\\%of the time, whereas the fully connected \(100%100\\%edge density\) problem requires around 774 ms and achieves a success rate of97%97\\%\.

![Refer to caption](https://arxiv.org/html/2607.02112v1/Figure4.png)Figure 4:\(a\) Time to solution and success rate for different edge density problems, with target Ising energy of90%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±1\\pm 1standard 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\.Edge DensityBest Ising EnergyHbSBBAWIM10%\\%\-803,331,77693\.10%93\.10\\%25%\\%\-1,201,604,48093\.59%93\.59\\%50%\\%\-1,491,890,81695\.52%95\.52\\%75%\\%\-1,535,851,13695\.74%95\.74\\%100%\\%\-1,339,657,85694\.52%94\.52\\%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\.

### 2\.4Number Partitioning

In this section, we focus on the number partitioning problem \(NPP\), which is one of Karp’s 21 NP\-complete problems\[Karp1972\]\. 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](https://arxiv.org/html/2607.02112#S2.F5)\(a\)\. The NPP has practical applications ranging from multiprocessor scheduling and minimizing the size and delay of the VLSI circuits\[coffman1991probabilistic,tsai1992asymptotic\], to public\-key cryptography\[merkle1978hiding\], and in choosing sides in a ball game\[hayes2002computing\]\. Furthermore, both the exact cover and the knapsack problems can be reduced linearly\[Karp1972,li2025efficient\]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​\(𝒜,ℬ\)=\|∑i∈𝒜ai−∑i∈ℬai\|E\(\\mathcal\{A\},\\mathcal\{B\}\)=\\left\|\\sum\_\{i\\in\\mathcal\{A\}\}a\_\{i\}\-\\sum\_\{i\\in\\mathcal\{B\}\}a\_\{i\}\\right\|\(4\)wherea1,a2,…\.aNa\_\{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∑ai\\sum a\_\{i\}is even\(odd\)\.

Lets=\+1s=\+1represent subset𝒜\\mathcal\{A\}ands=−1s=\-1represent subsetℬ\\mathcal\{B\}, then the above equation be written as

E​\(s\)=\|∑i=1Nai​si\|E\(s\)=\\left\|\\sum\_\{i=1\}^\{N\}a\_\{i\}s\_\{i\}\\right\|\(5\)So, the coupling terms for the Ising model for the NPP can be written asJi​j=ai​ajJ\_\{ij\}=a\_\{i\}a\_\{j\}with bias fieldh=0h=0\. This Hamiltonian for NPP is similar to the Mattis spin glass Hamiltonian\[mattis1976solvable,grass2016quantum\]\.

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 elementsJi​jJ\_\{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](https://arxiv.org/html/2607.02112#S2.F5)\(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 the99\.9%99\.9\\%approximate solution indicates that the difference is less than0\.1%0\.1\\%of the maximum possible value \(i\.e\., the total sum of all numbers\)\. In all cases, BAWIM outperforms the HbSB algorithm, achieving a100%100\\%success rate in most instances when considering the99\.9%99\.9\\%approximate solution\.

![Refer to caption](https://arxiv.org/html/2607.02112v1/Figure5.png)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\.5Sudoku

Sudoku is a famous puzzle game played on a9×99\\times 9grid, 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”\[resnick2014sudoku,mucke2024simple\]\. Beyond attracting puzzle enthusiasts, it has also become a frequent subject of study for mathematicians and computer scientists\[resnick2014sudoku,pignari2025efficient,mucke2024simple,shukla2025non,ercsey2012chaos,lynce2006sudoku\], likely owing to its simplicity and accessibility to a non\-technical audience\. For generalized grid sizes ofn2×n2n^\{2\}\\times n^\{2\}, the problem has been proven to be NP\-complete\[yato2003complexity\]\.

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 of81×9=72981\\times 9=729spins 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\[mediumSUDOKUSolvingWorlds\]\. Clues are embedded in the bias fieldhh, 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](https://arxiv.org/html/2607.02112#S2.F6)\(a\) illustrates the evolution of the Ising energy and the number of rule\-violating cells during the Sudoku\-solving process using BAWIM\. Figures[6](https://arxiv.org/html/2607.02112#S2.F6)\(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](https://arxiv.org/html/2607.02112#S2.F6)\(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\.

![Refer to caption](https://arxiv.org/html/2607.02112v1/Figure6.png)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\.6Estimation 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\.7Comparison 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\[inagaki2016coherent\]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\[honjo2021100\]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 approximately1×1061\\times 10^\{6\}carrier periods, while the more recent 100,000\-spin CIM used around 6,000 with a15%15\\%duty cycle\. In contrast, BAWIM operates with just 6 carrier periods per spin and with a50%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\[honjo2021100\]has a temperature coefficient of phase accumulation of1\.73×1071\.73\\times 10^\{7\}d​e​g/C∘deg/\{\{\}^\{\\circ\}\}C, meaning that a temperature change of just1∘​C1^\{\\circ\}Ccan 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±0\.05∘​C\\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 780d​e​g/∘C\{deg/^\{\\circ\}\}C, which is2\.21×1042\.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\.

## 3Conclusion

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\.

## 4Methods

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