A parallel magnetic tunnel junction-based probabilistic Ising processor for efficient quadratic optimization

Domain: biology
Published: 2026-03-30
Editorial outcome: minor_revision

Abstract
Solving computationally demanding combinatorial optimization problems using conventional computing architectures is slow and energy intensive. Quantum computing could improve optimization efficiency but remains at an early stage. Probabilistic computing offers a practical near-term approach to faster optimization through stochastic techniques.

Here, we experimentally demonstrate a scalable spin-transfer-torque-magnetic-tunnel-junction based probabilistic processor for efficiently solving all-to-all connected quadratic assignment problems. Our system integrates 144 compact spintronics tunable random number generators with a massively parallel architecture, achieving a high Monte-Carlo sampling throughput of 14.4 million flips per second. We co-design a parallel trial annealing scheme, and the integrated system achieves a 123× speedup with 98.3% energy savings over conventional Gibbs sampling, and a 3.2× speedup with 58.3% energy savings relative to the central processing unit implementation based on a compiled language. We further benchmark performance across graphics processing unit, and D-Wave quantum annealers, showing gains in solution quality, speed, and energy efficiency.

Introduction
Combinatorial optimization problems (COPs) arise in diverse domains1, including scheduling2, logistics3, circuit design4 and financial management5. Accelerating the solution of COPs can broadly impact numerous businesses. However, the end of Moore’s law has led to diminishing performance improvements in general-purpose processors6, prompting a shift toward the development of quantum computer7 and domain-specific probabilistic computing hardware8-10. While quantum computing has attracted attention as a potential revolution to optimization field11, current systems have yet to demonstrate practical speedups for industrial applications and remain constrained by complexity and high cost12. This has motivated a parallel exploration of classical physics based probabilistic computing13-15. These architectures increasingly leverage emerging technologies such as spintronics16-22, optics23, memristors24, field-programmable gate arrays (FPGAs)25-27, and other unconventional platforms28-30.

Given that massive random number generation is a key ingredient of probabilistic computing, magnetic tunnel junctions (MTJs) emerge as promising spintronic candidates due to their energy and area efficiency, rooted in their inherent stochastic switching behavior and tunable switching probabilities31-34. The statistical properties of the MTJ-generated random numbers have been characterized for security-oriented applications35. Spin-transfer torque magnetic tunnel junctions (STT-MTJs), the building blocks of magnetic random-access memory (MRAM), are well suited for constructing probabilistic computing hardware36. STT-MTJs offer distinct advantages such as superior device uniformity, high complementary metal-oxide-semiconductor (CMOS) compatibility, compact structure, controllable switching probability, and fast switching speed in the nanosecond regime. These attributes make STT-MTJ-based tunable random number generators (RNGs) strong candidates for large-scale hardware implementations of probabilistic computing systems37-39.

In this article, we present a parallel 144 STT-MTJ based probabilistic processor for solving challenging COPs on fully connected graphs. This processor leverages STT-MTJs for both true random number generation and Markov-Chain-Monte-Carlo (MCMC) sampling, reducing the latency, area, and energy associated with conventional random sources35, and enabling highly efficient intrinsic probabilistic sampling in physical devices. The quadratic assignment problem (QAP) is a fundamental NP-hard problem that models facility-location assignment tasks with facility pairwise interaction costs40. Due to its factorial search space and densely coupled structure, QAP is regarded as one of the most challenging problems in combinatorial optimization, with exact solvers scaling poorly beyond small instances41. In this work, we use QAP as a representative case study to investigate the potential computational advantages of probabilistic computing. A system implementing our MTJ-based parallel-trial architecture solves a 12-city QAP 123 times faster than conventional simulated annealing executed on the same hardware, while reducing energy consumption by 98.3%. Compared to a von Neumann type central processing unit (CPU, Intel i9-14900K) based implementation using a compiled language C++, the proposed parallel architecture achieves a 3.2× speedup and 58.3% energy savings. As the system scales, both performance and energy efficiency are expected to improve proportionally, underscoring the benefits of parallelism and the suitability for large-scale optimization. Beyond classical baselines, we also benchmark the solution quality and feasibility of our system on QAP instances against the state-of-the-art D-Wave quantum annealer with both Zephyr and Pegasus topologies.

Results
STT-MTJ based parallel probabilistic Ising machine

The solving strategy of probabilistic Ising machines (PIMs) involves mapping a COP to an Ising Hamiltonian, thereby connecting the optimal solution of the COP to the ground state of Ising energy landscape. In a lattice with N spins, where each spin is in states s = {s1,…,sN}, the Ising Hamiltonian (ℋ) is defined as ℋ(𝒔) = −∑N i<j Jij si sj −∑N i=1 hi si, where si and sj represent the states of the ith and jth spins, which can be either +1 or −1. Each ith spin is biased by a local field hi and interacts with the jth spin through a coupling Jij. The goal is to find the spin configuration that minimizes ℋ, i.e., the ground state. Unlike deterministic solvers, a PIM explores the solution space via stochastic spin updates, guided by a combination of energy minimization and sampling techniques such as Metropolis-Hastings and Gibbs sampling.

An STT-MTJ device consists of two ferromagnetic layers: a fixed reference layer with a stable magnetization, and a free layer whose magnetization can be switched by the STT effect49. These layers are separated by a thin magnesium oxide (MgO) insulating barrier, which enables quantum tunneling of spin polarized electrons50. When the magnetic moments of the two ferromagnetic layers are aligned parallel to each other, the device exhibits a low-resistance state, representing a spin state of +1; when they are antiparallel, it exhibits a high-resistance state, corresponding to a spin state of −1. With appropriately chosen pulse width and amplitude, STT induces stochastic switching, allowing the MTJ device to switch states with a controllable probability51.

Fig. 2a illustrates the system diagram of the STT-MTJ based PIM. The PIM comprises nine processing elements (PEs), each containing 16 STT-MTJ cells. Every 16 cells are driven by a 16-channel digital-to-analog converter (DAC) and sampled by a 16-channel analog-to-digital converter (ADC). Each STT-MTJ cell consists of one MTJ and one 3 kΩ sense resistor used for state readout. An FPGA interfaces with the STT-MTJ board to control the DACs and ADCs, and to compute the feedback signals for each MTJ cell at every annealing iteration. The feedback is based on the spin states read out by the ADCs and calculated based one energy minimization techniques, which typically involve a multiply-accumulate (MAC) operation between the coupling matrix J and the spin vector s. This architecture supports fully parallel operation of MTJ cells, providing flexibility to accommodate a wide range of algorithms that can accelerate from parallelism. A photograph of the STT-MTJ system is shown in Fig. 2b. The STT-MTJ devices are wire-bonded to red breakout boards, with each board hosting eight MTJs. Two breakout boards together with a DAC and an ADC form one PE. Custom FPGA Mezzanine Card (FMC) and FMC Plus (FMCP) to IPEX conversion boards are designed to enable high-speed data transmission between the FPGA and the MTJ board.

To achieve optimal problem-solving performance, a calibration process is performed after system setup. Based on the individual characteristics of each device, a linear transformation is applied to the input voltage amplitude of each MTJ. This transformation centralizes the probabilistic switching curves around zero and normalizes the slope at the origin to a value of 0.25. In practice, all MTJ cells are operated in parallel by fixing the DAC output pulse width to 4 µs, while sweeping the amplitude from 0.68 V to 1.34 V in steps of 4.5 mV. The choice of 4 µs perturb pulse width represents a practical trade-off between DAC/ADC speed and channel count in our proof-of-concept system. The resulting switching probability versus voltage curves for all MTJs are shown in Fig. 2c, where noticeable deviations and differences in slopes can be observed. The observed variation in MTJ switching probabilities mainly arises from fabrication-induced non-uniformities, such as variations in film thickness and etching precision. After characterizing the probabilistic switching behavior of each device, the curves are calibrated to align with a standard sigmoid function, as illustrated in Fig. 2d (see Supplementary Note 1 for details). Although the calibration cannot fully compensate for all device-to-device variations and non-idealities, it significantly improves performance, particularly by centralizing the switching curves around zero.

For large‐scale implementations (e.g., on the order of millions of MTJs), assigning distinct calibration parameters to every individual device would impose prohibitive hardware overhead. Importantly, wafer-level MRAM production has demonstrated excellent device-to-device uniformity in MTJ switching statistics. Based on this, it is reasonable to assume that only a small number of calibration parameter sets (e.g., one global set, or a few sets per process/wafer class) would be sufficient to capture the switching probability characteristics across all MTJs, without consuming significant memory resources. Furthermore, as the system scales, data transfer, particularly between memory and computing units, can substantially increase computation time. Benefiting from the CMOS compatibility, MTJ-based probabilistic computing systems can be integrated with emerging in-memory computing technologies such as MRAM52 and resistive random-access memory (RRAM)53 based platforms.

MTJ based parallel probabilistic Ising processor architecture

The digital annealer (DA) is a specialized hardware accelerator designed to solve combinatorial optimization problems by implementing a simulated annealing algorithm in a parallel-trial framework54,55. Each annealing iteration consists of two stages: a parallel-trial phase and a single-update phase. During the parallel-trial phase, the system evaluates N energy changes ΔE1, ΔE2…ΔEN in parallel, in analogy to quantum superposition, where ΔEi corresponds to the change in system energy if bipolar spin si is flipped to −si, while all other spins remain fixed. The acceptance for all proposed flips is then determined simultaneously using either the Metropolis-Hastings or Gibbs sampling criteria. In the subsequent single-update phase, the system randomly selects and commits one of the accepted spin-flip candidates, ensuring at most one spin changes its state per iteration to preserve the Markov chain dynamics.

In conventional CMOS based parallel-trial architecture, pseudo random number generators (PRNGs), such as linear feedback shift registers (LFSR) and lookup tables (LUT), are typically employed in the sampling process. For example, in the case of Gibbs sampling, the probability of accepting a spin flip follows a sigmoid-like function, 1/(1+exp(βΔE)), where β is the annealing parameter. This operation generally involves using a LUT to compute the sigmoid value and a 32-bit LFSR to generate a fix-point random number uniformly distributed between 0 and 1. A comparator is then used to determine the acceptance of the proposed flip by comparing the sigmoid result and the generated random number. PRNGs are also highly demanded in the single-update phase to randomly select one spin for flipping among all accepted candidates. Given the parallel nature of architecture and the large scale of COPs, area and energy consumption scaling can become significant.

In contrast, STT-MTJs inherently exhibit a sigmoid-like stochastic switching behavior that can be exploited for both the parallel-trial phase (MCMC sampling) and the single-update phase (random selection among accepted candidates). Rather than relying on bulky LUTs and LFSRs, sampling with MTJ is reduced to a simple operation: the computed value of βΔE is converted into voltage domain and applied to an STT-MTJ, whose probabilistic switching response yields the sampling result. Random number generation, required for spin selection in the update phase, is similarly achieved by biasing the MTJ at its 50-50 switching point, where thermal noise induces stochastic transitions with approximately 50% probability. It is estimated that an MTJ-based implementation can reduce transistor count by over three orders of magnitude (~1200×) and energy consumption by more than 70× compared to conventional CMOS based 32-bit LFSR + LUT implementation56. Moreover, unlike LFSRs57, the true random numbers generated by MTJs do not exhibit periodicity, regardless of runtime duration. The inherent and tunable randomness enhances the efficiency of the system to escape local minima and explore the energy landscape more effectively58.

Fig. 3 illustrates the MTJ-based parallel-trial accelerator, implemented as a tightly coupled FPGA-MTJ system. In each annealing iteration, the FPGA computes energy changes ΔE for all spin flip trials in parallel and converts them into MTJ voltage inputs through a linear transformation. These voltages are applied to MTJ cells, with each MTJ cell functioning as a local MCMC sampler. The analog outputs of MTJ cells are digitized by on-board ADCs and fed back to the FPGA. The FPGA then compares the digitized outputs against predefined thresholds to generate acceptance flags (denoted as flag in Fig. 3) for each trial. Subsequently, all MTJ cells are perturbed at their calibrated 50-50 switching point to generate true random bits, which are utilized by the update arbiter module.

The update arbiter module adopts a multi-layer structure to select a single spin among the pool of accepted candidates. Each input node in the arbiter carries acceptance flag and spin index. At each layer, a selector evaluates the two input nodes: if both have flags set to 1, one spin index is randomly forwarded based on MTJ-generated random bits; if only one flag is 1, that spin index is passed through; and if both flags are 0, a null index with flag 0 propagates. This arbitration process continues across layers until a single spin index is selected for update, ensuring that at most one spin is flipped per annealing iteration.

Evaluation on quadratic assignment problems

We evaluate the performance of our STT-MTJ based parallel-trial accelerator using problem instances from the tinyQAP dataset59. This dataset includes QAP instances with sizes ranging from 3 to 12 facilities and locations (see Supplementary Note 2 for additional tests on large instances). The mapping from the QAP to the Ising model follows the quadratic scaling rule: each QAP instance with n facilities and locations is transformed into an Ising model comprising n2 spins. As a result, the corresponding Ising graphs span from fully connected 9 to 144 spins.

In the conventional simulated annealing (SA) scheme, only one spin-flip proposal is evaluated and updated in an annealing iteration. Therefore, using a single MTJ to sequentially update all spins is equivalent to using n2 MTJs to update all spins, without any loss in speed or energy consumption. Based on this, we use a single STT-MTJ for spin update in our system for SA. However, in the parallel-trial annealing (PA) scheme, where all spin-flip proposals must be evaluated simultaneously. The number of MTJ devices should match the number of spins in the Ising model to fully exploit parallelism. Leveraging the massive parallel architecture of our STT-MTJ based system, we utilize 2n2−1 MTJs in parallel for spin update in the PA scheme, where n2 MTJs are used for MCMC sampling, and n2−1 MTJs function as random number generators for the arbiter module during the single-update phase.

We benchmark our MTJ-based probabilistic Ising processor against digital platforms, including a high-performance Intel i9-14900K CPU and an NVIDIA H100 GPU, on QAP instances. On the CPU, we implement simulated annealing in both an interpreted environment (MATLAB) and a compiled environment (C++) to provide a fair comparison across different software execution paradigms. On the GPU, we implemented replicated simulated annealing (replicated SA) in C++ environment. Specifically, we implement the ΔE update computations within CUDA threads and assign 1024 independent SA replicas to separate CUDA thread blocks. Further GPU implementation details are provided in Supplementary Note 3.

Fig. 4a presents the median time to reach 95% of the optimal solution (TTS95) as a function of problem scale, measured by the number of spins in the Ising model. All computing platforms are evaluated under an identical annealing schedule to ensure a controlled comparison. Relative to MTJ-SA, the MTJ-PA scheme exhibits a substantial acceleration, achieving speedup factors of up to 123×. When benchmarked against CPU-based SA, MTJ-PA provides improvements of up to 6.6× over the interpreted-language MATLAB implementation and 3.2× over the compiled-language C++ implementation. In comparison with GPU-based replicated SA, MTJ-PA attains up to 8.9× acceleration. It should be noted that, in our current implementation, the GPU and CPU exhibit comparable performance. This behavior arises because the problem sizes used in our present evaluation are not sufficiently large for the massive parallelism of GPU to yield a significant advantage. For relatively small instances, the overhead associated with GPU kernel launches and data movement offsets the benefits of parallel execution, leading to performance similar to that of the CPU.

When comparing MTJ-SA with CPU-SA, we observe that the MTJ-based SA operates more slowly than its CPU counterpart, primarily due to peripheral circuit limitations in our proof-of-concept system, most notably the constrained DAC/ADC throughput (see Supplementary Note 4 for a detailed system-latency analysis). Despite this limitation, the MTJ-based PA still substantially accelerates the annealing process by leveraging both the algorithmic advantage of a higher spin-flip acceptance rate and the parallel architecture enabled by tunable MTJ-based random number generators. To further enhance the performance of MTJ-based probabilistic computing systems, future implementations can leverage tighter MTJ-CMOS integration through back-end-of-line (BEOL) technologies to achieve more compact architecture. By co-integrating MTJs with on-chip high-speed DACs, sense amplifiers, and local control logic, the overall system latency can be significantly reduced. With MTJs operating at their intrinsic switching limit, on the order of 250 million flips per second per device38,60, the achievable acceleration could approach four orders of magnitude compared to CPU-based implementations. Such a speedup is also expected to improve solution quality within the given runtime budget. It is also worth noting that we use a simple linearly increasing inverse temperature schedule, without any fine-tuning or specialized optimization. Incorporating more advanced algorithmic techniques61 into our architecture could further improve performance.

Fig. 4b presents the median energy-to-solution (ETS) scaling, including energy consumption from FPGA, DACs, ADCs and MTJs as a function of QAP problem size (represented by the number of spins). The energy-to-solution is estimated by multiplying the system power (see Method for details) by the corresponding time-to-solution. Although MTJ-PA consumes more instantaneous power than MTJ-SA due to its parallel architecture (see Supplementary Note 5 for details), it still achieves lower overall energy consumption because of the significant reduction in TTS. The energy saving factor is defined as 1−(ETS_MTJ-PA /ETS_other scheme), where ETS_MTJ-PA and ETS_other scheme represent the energy-to-solution of MTJ-PA and other schemes. The energy-saving factor between MTJ-PA and MTJ-SA increases with the level of parallelism, reaching 98.3% in our system. When compared to CPU-SA, the MTJ-PA achieves 87.1% energy savings over the interpreted-language MATLAB implementation and 58.3% over the compiled-language C++ implementation, demonstrating the significant efficiency advantage of MTJ-based stochastic computing.

Benchmark digital, probabilistic and quantum computing

To assess the performance of conventional digital computing, classical physics-based probabilistic computing, and quantum computing in solving complex optimization problems, we benchmarked representative platforms on a series of fully connected QAPs with sizes ranging from 3 to 12. For digital computing, we implement a simulated annealing algorithm on a CPU with 5000 sweeps, where each sweep consists of attempted updates for every spins. For probabilistic computing, we employ a parallel-trial annealing scheme on our STT-MTJ-based probabilistic Ising machine, running for 5000 iterations, with each iteration updating at most one spin. For quantum computing, we benchmark with the state-of-the-art D-Wave quantum annealer.

The D-Wave annealer employs a process known as quantum annealing to solve optimization problems by finding the lowest-energy configuration of a system of superconducting qubits. Each optimization problem is first formulated as an energy landscape, defined by setting the qubit biases and coupler strengths that represent the objective function of the problem. The system begins with all qubits in a quantum superposition of possible states. During the annealing process, the Hamiltonian parameters are gradually evolved so that quantum fluctuations and entanglement guide the system through many possible configurations. As the anneal progresses, the influence of the driver Hamiltonian is reduced while that of the problem Hamiltonian is increased, allowing the system to evolve naturally toward low-energy states62. Because the connectivity of the D-Wave hardware graph (e.g., Pegasus and Zephyr topologies) does not always match the connectivity of the problem graph, minor embedding is required to map logical variables and couplings onto the available physical qubits. In this process, a single logical variable may be represented by a chain of physically connected qubits, with the chain strength tuned to preserve variable consistency during annealing. At the end of the anneal, typically within time scale of microseconds, the quantum superposition collapses, and the qubits settle into a single classical configuration corresponding to the minimum (or near-minimum) energy of the encoded problem, representing the final solution. As a commercially available quantum computing platform, D-Wave has demonstrated the ability to address certain real-world optimization problems. Thus, in this work, we use D-Wave as a representative quantum annealer to provide a hardware benchmark alongside the digital CPU and our MTJ-based probabilistic annealer. Specifically, we use the Advantage_system6.4 solver with Pegasus topology63 and the Advantage2_prototype1.6 solver with Zephyr topology64 (Detailed implementation parameters are included in Supplementary Note 6). A long annealing time of 2000 µs is configured for both quantum solvers to maximize solution quality for the optimization problems (see Supplementary Note 7 for a systematic analysis of annealing time).

Fig. 5a presents the scaling of solution quality as a function of QAP problem size. All solutions are normalized with respect to the best-known solutions obtained using the commercial solver Gurobi (see Supplementary Note 8 for Gurobi implantation details). Only feasible solutions, those that satisfy all QAP constraints, are included in this analysis. As the problem size increases, MTJ-PA consistently delivers high quality solutions closely track to the optimal baseline. CPU-SA performs slightly worse than MTJ-PA but remains comparably close to the optimum. In the quantum computing regime, both D-Wave Zephyr and D-Wave Pegasus can find feasible solutions for QAP instances with fewer than four facilities and locations. Zephyr slightly outperforms Pegasus in terms of average solution quality but both solvers exhibit substantial variability, as indicated by large standard deviations. Beyond the scale of 5, both the Pegasus and Zephyr solvers fail to return a single feasible solution across 100 runs. As shown in Fig. 5b, the probability of finding feasible solutions drops sharply for both quantum solvers and becomes effectively zero from size 5 onward. In contrast, both MTJ-PA and CPU-SA consistently yield feasible solutions across the entire dataset. Although the performance of quantum solvers may be improved with more advanced embedding and mapping strategies65, they remain at an early stage of development and are still far from demonstrating practical advantages for industry-scale applications. In contrast, STT-MTJ-based probabilistic computing offers a promising and scalable approach for accelerating computationally hard optimization problems.

Methods
Thin film growth and device fabrication. Thin film stacks with the structure: substrate/[W (3)/Ru (10)]₂/W (3)/Pt (3)/Co (0.25)/Pt (0.2)/[Co (0.25)/Pt (0.5)]5/Co (0.6)/Ru (0.85)/Co (0.6)/Pt (0.2)/Co (0.3)/Pt (0.2)/Co (0.5)/W (0.3)/CoFeB (0.9)/MgO (0.85)/CoFeB (1.2)/W (0.4)/CoFeB (0.8)/MgO (0.7)/Ta (3)/Ru (7)/Ta (5) (thicknesses in nanometers) were deposited on Si substrates with a 300 nm thermal oxide layer. The deposition was carried out at room temperature using direct current (DC) magnetron sputtering for metallic layers and radio frequency (RF) magnetron sputtering for MgO, under a base pressure below 2 × 10⁻8 Torr. To fabricate the MTJ devices, 10 µm-wide bottom electrode structures were first patterned via photolithography and Ar ion milling. MTJ pillars with ~90nm diameter were then defined using electron beam lithography. Following ion milling, an Si3N4 encapsulation layer was deposited in situ using RF magnetron sputtering, without breaking vacuum. Finally, 10 µm-wide top electrode structures were patterned using photolithography, and the top electrode layers, Ta (5 nm)/Cu (40 nm), were deposited using DC magnetron sputtering.

MTJ carrier printed circuit board (PCB) design. Sixteen 16-channel DACs (AD5767) are used to generate analog input to each MTJ, and sixteen 16-channel ADCs (MAX11131) are used to sample the output of each MTJ cell and send to the FPGA (AMD Virtex UltraScale+ FPGA VCU118 Evaluation Kit) for feedback signal calculation. The polarity of Vin dictates the current direction such that the reset and perturb processes can be realized. Both DACs and ADCs interface with the FPGA through serial peripheral interface (SPI) protocol via IPEX hardware connection. The RTL logic in FPGA is programmed and implemented in Vivado 2023.2 platform using Verilog. In our setup, the main clock frequency in the FPGA is set to 100 MHz. For QAP instances with scale > 8, we use two MTJ carrier PCBs connected to the FPGA via IPEX and controlled in parallel by the VCU118. One PCB serves solely as a random number generator bank, with its MTJs biased at the 50-50 switching probability point to supply random numbers to the update-arbiter module.

Matrix-vector multiplication in FPGA. In the parallel-trial phase, each trial assumes a single spin flip while all other spins remain fixed. Under this assumption, the energy change term ΔEi associated with spin si flip to −si is given by: ΔEi = −2siLi, where Li(𝒔) = hi +∑N j=1 Jij sj. All the initial Li values are precalculated and stored in registers, as illustrated in Fig. 3 blue block. When a spin state is updated, each Li is incrementally updated by adding its corresponding change Li. For example, when sj is flipped to −sj, the resulting increment Li = Jij (−2sj).

TTS estimation. In the FPGA design, an integrated logic analyzer (ILA) is implemented to probe the current state of the finite state machine (FSM) that controls the annealing process. The duration of a single annealing iteration can be measured using the ILA, and multiplying this duration by the number of iterations required to reach a solution yields the total TTS.

System power estimation. For the MTJ-based system, the total power consumption was estimated as the sum of the FPGA, ADC/DAC, and MTJ contributions. The FPGA power was obtained from the Vivado Power Report generated after design implementation. The ADC and DAC power values were calculated from their respective supply voltage and current specifications. The MTJ power was estimated by summing the consumption during the reset and perturb phases. For the CPU-based implementation, the power consumption is estimated by multiplying the CPU utilization rate of the program (MATLAB or C++ executable) during runtime by the nominal power rating of the processor. For the GPU-based implementation, we estimate a light-load power baseline of 90 W.

Data availability. The data that support the findings of this study are available within the article and the Supplementary Information.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript proposes an innovative hardware–algorithm co-design approach for solving quadratic assignment problems, which is both creative and potentially highly impactful. I recommend it for publication. By integrating STT-MTJs as true random number generators with a parallel-trial annealing algorithm, the authors have greatly reduced runtime and energy consumption. While this work achieves significant results, I recommend the following revisions to strengthen its rigor, clarity, and impact:\n(1) I understand the acceleration benefit of MTJs in generating single random numbers and in p-bit operations. However, in Fig. 4a, the SA running on a CPU is faster than the SA running on the Ising machine built in this work. Is this due to FPGA signal transmission latency or MTJ size limitations? Please provide further analysis. The manuscript also states that MTJ-based PA runs faster than CPU-based SA — is this advantage primarily due to the algorithm or the MTJs? Can the authors break down the contributions from each? Based on current results, PA on this Ising machine likely takes longer than PA on a CPU. In the future, how should MTJs be integrated with FPGA or ASIC chips to surpass CPU runtime performance?\n(2) In the parallel-trial annealing algorithm described in the manuscript, the parallel-trial phase involves applying voltages to the MTJs for sampling to identify candidate nodes, followed by randomly selecting one candidate to transfer to the FPGA for the actual update. This represents a compute–memory separation strategy that does not exploit the non-volatility of the MTJs. In this architecture, could the MTJs be replaced by other volatile p-bit devices? If so, would different types of volatile p-bits exhibit differences in sampling speed, throughput, and noise characteristics, and how would these differences affect the final solution quality? If non-volatility is not utilized, do MTJs still offer unique advantages—such as more favorable sigmoid probability characteristics, better thermal stability, or CMOS compatibility—over other volatile p-bit implementations?\n(3) The manuscript emphasizes high sampling throughput (14.4M flips/s) but does not analyze system latency, such as DAC/ADC conversion delays, FPGA computation delays (especially MAC operations), or data communication delays. Were these latencies considered when calculating throughput? How was 14.4M flips/s derived — from a theoretical value (single MTJ flip rate × number of devices) or from measured values including read/write and processing delays? 14.4M flips/s corresponds to 10 µs per flip per device. Why was this operating time chosen instead of approaching the device’s physical limit (as stated in the manuscript, 250M flips/s, 4 ns/flip)? Has 14.4M flips/s reached or exceeded the throughput of current CPUs for equivalent MCMC sampling tasks? If so, have the authors considered using larger MTJs (trading density and speed for simpler fabrication and higher reliability) for certain practical applications?\n(4) The results show that the proposed STT-MTJ parallel probabilistic Ising machine can stably achieve solution quality close to the optimum (~95%) for QAP instances. Within the current hardware and algorithm framework, how could solution quality be improved from 95% to the exact optimum? For example, by extending annealing time, increasing iteration count, introducing multi-start parallel search, increasing the number of MTJs, or reducing system latency (to increase flip frequency)? If this is not achievable with the current framework, the authors should discuss possible strategies for progressing from near-optimal to exact solutions.\n(5) The manuscript reports a noticeable variation in MTJ switching probabilities (Fig. 2c) and compensates for it through calibration — a valuable approach to addressing device uniformity issues. From a device perspective, what might be the underlying causes of such variation? Given the relatively narrow voltage tuning range and device-to-device inconsistency, could operational stability be affected (e.g., probability curve shifts over time or temperature, while stored σ and μ remain unchanged)? Could changes in ambient magnetic fields impact global applicability? Moreover, the implications for large-scale scaling (e.g., millions of devices) are not discussed. Would storing σ and μ for all devices require significant FPGA memory, and would large-scale data transfer substantially increase computation time?\n(6) In Supplementary Tables S1 and S2, the authors present power consumption values for FPGA, ADC+DAC, and MTJ components under different system scales. Were these values measured experimentally or calculated from device parameters and duty cycles? Providing more detail on the measurement or estimation methodology would enhance clarity.\n(7) On page 5, second paragraph, the full form of the abbreviation “PIM” is not given.", "reviewer_2": "In this work, the authors present an experimental demonstration of spin-transfer-torque magnetic tunnel junction (STT-MTJ) based probabilistic processor for solving quadratic assignment problems (QAPs). The integration of 144 MTJs with an FPGA represents a significant step forward in the field of spintronics-based probabilistic computing. The work is novel, well-executed, and the results are promising. However, several concerns must be addressed to strengthen the manuscript's clarity before it can be accepted for publication.\n(1) \"PIM\" (Probabilistic Ising Machine) is introduced without prior definition. Please introduce and define this term clearly for readers who may not be specialists in this specific sub-field.\n(2) Regarding Fig. 1d, the TSP energy landscape is depicted as smooth, while QAP is rugged. Is this true for all TSP instances, or just a selected case? The authors should clarify this point and provide appropriate citations to literature that discusses the comparative complexity of TSP vs. QAP landscapes to avoid over-generalization.\n(3) Regarding the description of emerging spintronic technologies for probabilistic computing, it would be beneficial to include references on related works, such as: Nano Lett. 2024, 24, 5420−5428 (https://doi.org/10.1021/acs.nanolett.3c04820).\n(4) The experiments are limited to 12-city QAPs, which may be too small to convincingly demonstrate scalability. What is the maximum scale (number of cities)? At what number of cities does this proposed system fail to achieve the optimal solution? Please test larger instances and report when the performance degrades or fails to find the optimum.\n(5) Regarding Fig. 4c, the comparison is conducted between MTJ-PA and CPU-SA. Why was there no comparison with CPU-PA? To demonstrate the superiority of the MTJ-integrated system, a comparison between MTJ-PA and CPU-PA is essential.\n(6) The authors have mentioned “compared to CPU-SA, the MTJ-PA configuration achieves energy savings of up to 82%”. How was the energy consumption of CPU-SA measured? This manuscript does not mention how the energy/power consumption of CPU-SA was calculated/measured. Please provide a detailed methodology.", "reviewer_4": "Classical probabilistic computing, in particular for topics of monte carlo sampling, is certainly of interest in terms of computational technologies. The comparison to simulated annealing is good to show. The comparison against D-Wave quantum annealers is a somewhat strange comparison choice in my view -- typically, especially for these types of optimization problems, simulated annealing performs much better than D-Wave quantum annealers.\nMost importantly, having the ground state optimal solution found by Gurobi is an excellent part of the study -- this is crucial to show if the evaluated methods can reach the ground-state energy.\nSpecific comments:\n1. The methods used in the D-Wave QPU runs needs to be specified in complete detail. Namely, what minor embeddings were used, chain strength, broken chain repair algorithm. Also, importantly the total number of anneals measured needs to be stated (and whether the solution quality is an average, or the best). Also, the total number of logical variables in the original QUBO needs to be stated for complete methods.\n2. Figure 1d: it should be clarified that this is purely a schematic, not actual data, it is somewhat confusing otherwise. Line 93 on page 5 states that \"By contrast, the TSP plot looks much more like a smooth bowl \", which is confusing to me -- is this an empirical statement? I think it is somewhat hard to correctly interpret the search space of a high dimensional discrete combinatorial optimization search space. This portion of the manuscript needs to be clarified as to whether this is actual data, or is a conceptual interpretation of the search spaces.\n3. I would generally say that the methods description, and properties, of the D-Wave quantum annealing devices is not very sufficient -- it would be good to cite and describe the technology ideas, how it works, why it may be reasonable to use in this context, etc. Maybe also why minor embedding is required, maybe also how much of the QPU chip was used in the embedding process.\n4. The CPU compute time used by Gurobi should be reported as well.\n5. One overall comment I have on the aims and claims of this study is that in order to really make a state of the art methods comparison between a hardware probabilistic sampler and classical digital computers would be to also run a GPU comparison -- highly optimized QPU simulated annealing implementations can be much more efficient than CPU implementations.\n6. Line 243: \"For digital computing, we implement a simulated annealing algorithm on a CPU with 5000 sweeps\" For this type of computing paradigm comparison, the exact implementation really matters. Therefore, in this case, I think the authors should give more details on the programming language used, ideally a specific library that was used. Ideally this would be a compiled language, something along the lines of C, Julia, Rust, in order for the comparison to be very rigorous.\n7. Line 173: in regards to the random number generation, this point is interesting because thermodynamic based random number generation is interesting. However, it is not clear in this study what exactly is the advantage gained by this particular random number generator. Generally, I would argue that a dedicated analysis of the random number generator would be good (e.g. min entropy, and general random number generator tests such as from Diehard and the NIST testsuite) -- and then, showing what exactly is \"high-quality\" (to quote the manuscript) about this randomness.\n8. Selection of the annealing time on the D-Wave QPUs: a choice if 2000 microseconds is reasonable, generally it should result in the lowest sampled energy. However, it is possible for some types of optimization problems, a slightly shorter annealing time is better. Ideally, one would perform some sort of systematic annealing time analysis. But, I do not think this crucial for the article to be published.\n9. Line 80, page 4: \"function becomes equivalent to finding the ground state spin configuration\": is there any potential degeneracy in this formulation? This sentence implies that there is only ever one unique solution.\nOverall I would not recommend this manuscript for publication in its current form -- it requires revisions to clarify the methods, and some of the plots. In terms of general community interest, I would estimate that this study is not of considerable interest to the quantum computing community. But, I do think the most interesting part of the study is the hardware description of the probabilistic Ising model sampler."}, "author_response": "Response to the First Reviewer\nReviewer #1 (Remarks to the Author):\nThe manuscript proposes an innovative hardware–algorithm co-design approach for solving quadratic assignment problems, which is both creative and potentially highly impactful. I recommend it for publication. By integrating STT-MTJs as true random number generators with a parallel-trial annealing algorithm, the authors have greatly reduced runtime and energy consumption. While this work achieves significant results, I recommend the following revisions to strengthen its rigor, clarity, and impact:\nAuthors’ response: We sincerely thank the reviewer for the positive evaluation. Please find our detailed, point-by-point responses below.\n(1) I understand the acceleration benefit of MTJs in generating single random numbers and in p-bit operations. However, in Fig. 4a, the SA running on a CPU is faster than the SA running on the Ising machine built in this work. Is this due to FPGA signal transmission latency or MTJ size limitations? Please provide further analysis.\nAuthors’ response: We thank the reviewer for the constructive suggestion. The speed gap in simulated annealing between the CPU and our MTJ-based system is primarily limited by the DAC/ADC operation speed and FPGA signal transmission latency, rather than by the MTJ size. There are three dominant factors:\n1. DAC slew rate. To balance the trade-off between the number of channels and output speed, we used the LTC2688 (Analog Devices), a 16-channel precision DAC with a slew rate of ~1 V/µs. To ensure stable pulse generation and avoid transient voltage, the pulse width in our experiments was set to approximately 4 µs, rather than the intrinsic MTJ switching limit of ~4 ns (as shown in Fig. P2). This conservative setting ensures accurate biasing but limits the spin update rate.\n2. ADC sampling rate. The analog readout after each perturb pulse is digitized by multi-channel ADCs (MAX11131), whose sampling and conversion times are in the order of microseconds per channel.\n3. Signal-transmission latency between FPGA and DAC/ADC interface. The data exchange between the FPGA and the DAC/ADC boards is handled through a serial peripheral interface (SPI). Each spin-update cycle requires serial command and data transmission for both bias writing and readout retrieval. The limited SPI clock frequency (50 MHz) and long serial command words constrain the update speed.\nTo break down the time consumption, each MTJ based spin update comprises 5 stages: multiply accumulation (MAC), computation of the input voltage for each MTJ cell (CAL_Vin), DAC operation, ADC sampling, and decode ADC output (Decode). The breakdown in Fig. P1 illustrates the time consumption of each process in the simulated annealing operation on our MTJ-based system, revealing that the DAC and ADC operations dominate the total runtime. The duration of each stage is measured using the in-built Integrated Logic Analyzer (ILA) in the FPGA.\nFig. P1| Time analysis of MTJ system in simulated annealing on TinyQAP12\nWe further verify that the operation speed is not limited by the intrinsic MTJ device physics such as MTJ pillar size, through short pulse switching measurements. As shown in Fig. P2, the full probabilistic switching curve can be obtained using a pulse width as short as 4 ns, demonstrating that the MTJ devices can respond within nanosecond timescales. This result confirms that the device switching speed is not the bottleneck in the current system; rather, the overall latency arises from peripheral circuit and communication constraints.\nWe have added a concise explanation to the Main Text (page 11) and a detailed system latency analysis to Supplementary Note 4. The explanation in the Main Text reads “When comparing MTJ-SA with CPU-SA, we observe that the MTJ-based SA operates more slowly than its CPU counterpart, primarily due to peripheral circuit limitations in our proof-of-concept system, most notably the constrained DAC/ADC throughput (see Supplementary Note 4 for a detailed system-latency analysis).”\nThe manuscript also states that MTJ-based PA runs faster than CPU-based SA — is this advantage primarily due to the algorithm or the MTJs? Can the authors break down the contributions from each?\nAuthors’ response: The advantage of MTJ-based parallel-trial annealing (PA) over CPU-based simulated annealing (SA) arises from the algorithm, the MTJ device, and the parallel architecture.\n1. Algorithmic contribution—increase spin flip acceptance probability: Unlike simulated annealing (SA), which sequentially evaluates spin states and accepts updates one at a time, parallel-trial annealing (PA) identifies all spins that satisfy the update condition in parallel and selects among them. This parallel evaluation eliminates the idle (non-updating) iterations inherent in conventional SA, thereby substantially reducing redundant computations.\n2. MTJ device contribution—massively parallel, tunable random number generators. MTJ-based tunable true random number generators offer small footprints, good uniformity, and high energy efficiency, enabling large-scale implementations. This supports parallel random number generation for each spin flip trial evaluation block.\n3. Parallel architecture contribution—facilitate parallel evaluation of the spin flip trials.\nTogether, these three factors allow MTJ-based PA to outperform CPU-based SA. We have summarized this analysis to the Main text (page 11) as: “Despite this limitation, the MTJ-based PA still substantially accelerates the annealing process by leveraging both the algorithmic advantage of a higher spin-flip acceptance rate and the parallel architecture enabled by tunable MTJ-based random number generators.”\nBased on current results, PA on this Ising machine likely takes longer than PA on a CPU. In the future, how should MTJs be integrated with FPGA or ASIC chips to surpass CPU runtime performance?\nAuthors’ response: In a single-threaded CPU implementation, the parallel-trial annealing algorithm can only be emulated sequentially, as the CPU executes one instruction at a time. As a result, for a problem with n spins, the CPU based PA requires O(n) time per iteration to evaluate all spin flip proposals, which is effectively equivalent to one simulated annealing (SA) sweep. In contrast, on parallel hardware, all spins can be evaluated simultaneously, resulting in a time complexity of O(1) per iteration.\nAs shown in Fig. P3, we compare the time-to-solution of CPU-based SA, CPU-based PA, and MTJ-based SA and PA. Overall, the CPU-based PA exhibits a performance comparable to that of the CPU-based SA. For small problem sizes (e.g., scale below 7), the CPU-SA performs slightly better compared to CPU-PA because the CPU-PA implementation includes an additional arbiter process in each iteration to stochastically select among update candidates, which slightly increases the time consumption per iteration. However, starting from problem scale 8, CPU-PA begins to show modest improvement due to its ability to increase spin flip proposal acceptance probability. MTJ-based PA continues to outperform CPU-PA as the problem scales up, benefiting from the hardware-level parallelism inherent to the architecture, where all spin updates are evaluated concurrently.\nIn our prototype, MTJ-PA already outperforms CPU baseline on the quadratic assignment problem benchmark. To further enhance the performance of MTJ-based probabilistic computing systems, future implementations can leverage tighter MTJ-CMOS integration through back-end-of-line (BEOL) technologies to achieve a more compact architecture. By co-integrating MTJs with on-chip high-speed DACs, sense amplifiers, and local control logic, the overall system latency can be significantly reduced, enabling faster and more energy-efficient operation. We have revised and clarified the discussion on potential future performance enhancements in the Main Text (page 12).\nFig. P3| Time-to-solution comparison among CPU based PA, SA and MTJ based PA, SA.\n(2) In the parallel-trial annealing algorithm described in the manuscript, the parallel-trial phase involves applying voltages to the MTJs for sampling to identify candidate nodes, followed by randomly selecting one candidate to transfer to the FPGA for the actual update. This represents a compute–memory separation strategy that does not exploit the non-volatility of the MTJs. In this architecture, could the MTJs be replaced by other volatile p-bit devices? If so, would different types of volatile p-bits exhibit differences in sampling speed, throughput, and noise characteristics, and how would these differences affect the final solution quality? If non-volatility is not utilized, do MTJs still offer unique advantages—such as more favorable sigmoid probability characteristics, better thermal stability, or CMOS compatibility—over other volatile p-bit implementations?\nAuthors’ response: Thanks for such an intuitive question. In the current architecture, the non-volatility of MTJs is not explicitly utilized during the parallel-trial phase. In principle, they could be replaced by other volatile p-bit devices, such as volatile memristors [Woo, K.S. et al., Nat. Commun. 13, 5762 (2022)] [Rhee, H et al., Nat. Commun. 14, 7199 (2023)]. Depending on specific device, p-bits can exhibit differences in sampling speed, throughput and noise characteristics.\nSampling speed and throughput determine how many iterations can be executed within a given time window, thereby directly affecting the achievable solution quality. Device switching response sensitivity to input electrical pulses and noise characteristics jointly shape the effective sigmoid probability curve, influencing both its precision and stability. Devices exhibiting more controllable sigmoid transfer functions, that is, more finely resolved switching probability curves, can provide more accurate sampling of the Boltzmann distribution, improving convergence. In contrast, excessive temporal noise, device-to-device variability, or non-ideal sigmoid characteristics may distort the probabilistic sampling process and ultimately degrade the convergence quality of the solution.\nEven though non-volatility is not directly exploited in the present system, MTJs still offer several distinct advantages over other volatile p-bit implementations. These include\n1. Superior device uniformity, which are critical for large-scale parallel systems requiring millions of stochastic elements.\n2. High CMOS compatibility, enabling straightforward scaling and dense interconnection with peripheral control circuitry.\n3. Controllable switching probability and fast switching speed in the nanosecond regime, allowing precise tuning of probabilistic behavior and high-throughput sampling.\nTherefore, while other volatile p-bit devices could in principle perform similar probabilistic sampling functions, MTJs remain a robust and manufacturable candidate for scalable hardware implementation of probabilistic computing systems. The above discussion has been incorporated into the introduction section to highlight the distinct advantages of MTJ-based probabilistic computing systems.\n(3) The manuscript emphasizes high sampling throughput (14.4M flips/s) but does not analyze system latency, such as DAC/ADC conversion delays, FPGA computation delays (especially MAC operations), or data communication delays. Were these latencies considered when calculating throughput? How was 14.4M flips/s derived — from a theoretical value (single MTJ flip rate × number of devices) or from measured values including read/write and processing delays? 14.4M flips/s corresponds to 10 µs per flip per device.\nAuthors’ response: The estimated throughput of 14.4 M flips/s includes all system-level latencies. Specifically, it was measured implementing an Integrated Logic Analyzer (ILA) embedded in the FPGA.\nIn our implementation, the entire system operation is governed by a finite state machine (FSM), which sequentially controls each stage of the annealing process, including multiply–accumulate (MAC) operations, MTJ input-voltage calculation (CAL_Vin), DAC operation, ADC sampling, and spin decoding. The ILA captures the “r_current_state” signal in real time (measured in cycles) as the FSM transitions through each phase, enabling a direct measurement of the total cycle count for one complete operation.\nBy multiplying the measured cycle count by the clock period (100 MHz, 10 ns per cycle), the corresponding real-time duration can be calculated. Since the FSM inherently covers all intermediate steps and waiting periods, the throughput value derived from the ILA measurement fully accounts for DAC/ADC conversion delays, data communication overhead, and FPGA computation time. As every MTJ follows same FSM sequence, the total sampling rate is calculated as N /T , where N is the number of MTJ MTJ cycle MTJ cells operating in parallel, and T is the measured total cycle duration per FSM loop. A representative cycle latency breakdown for a 12-scale QAP instance is summarized in Table P1, showing the time consumption (in clock cycles) of each Markov-Chain-Monte-Carlo (MCMC) sampling stage.\nIn response to the reviewer’s comment, we have included the detailed discussion on system latency to Supplementary Note 4.\nTable P1 | Time consumption of each MCMC sampling stage.\nr_current_state Time (cycles) Time (µs)\nMAC 153 1.53\nCAL_Vin 9 0.09\nDAC operation 772 7.72\nADC sampling 62 0.62\nSpin decode 6 0.06\nTotal 1002 10.02\nWhy was this operating time chosen instead of approaching the device’s physical limit (as stated in the manuscript, 250M flips/s, 4 ns/flip)? Has 14.4M flips/s reached or exceeded the throughput of current CPUs for equivalent MCMC sampling tasks? If so, have the authors considered using larger MTJs (trading density and speed for simpler fabrication and higher reliability) for certain practical applications?\nAuthors’ response: The choice of 10 µs per flip instead of operating at the device’s physical limit of 4 ns per flip reflects a practical trade-off between DAC/ADC speed and channel count in our proof-of-concept system. In general, high-speed DACs and ADCs offer limited channel numbers, whereas multi-channel devices provide higher parallelism but operate at lower update rates. To support simultaneous biasing and readout of multiple MTJ cells, we selected multi-channel DACs/ADCs (LTC2688 and MAX11131), which prioritize channel density and measurement stability over speed. The 4 ns/flip, as shown in Fig. P2 represents the intrinsic MTJ switching limit, characterized independently using a Keysight 81134A pulse generator, confirming that the device-level stochastic response can operate at hundreds of megahertz (≈250 M flips/s). However, this speed cannot be fully utilized in the current setup due to the limitations of the peripheral circuits.\nCompared to modern CPU implementations, such as that reported by Aadit et al. [Nat. Electron. 5, 460–468 (2022)], which demonstrate a sampling throughput of approximately 0.2 M flips/s, our proof-of-concept MTJ-based system achieves a higher throughput of 14.4 M flips/s by leveraging hardware-level parallelism. We agree with the reviewer that, in future implementations, when MTJs are co-integrated with high-speed DACs and on-chip sense amplifiers, the fabrication process should be further optimized to achieve an optimal trade-off between the intrinsic MTJ switching speed, the operating speed of peripheral circuits, and device uniformity for practical applications.\nIn response, we added a brief explanation on our choice of 4 µs perturb pulse width to the Main text (page 7) as “The choice of 4 µs perturb pulse width represents a practical trade-off between DAC/ADC speed and channel count in our proof-of-concept system.”\n(4) The results show that the proposed STT-MTJ parallel probabilistic Ising machine can stably achieve solution quality close to the optimum (~95%) for QAP instances. Within the current hardware and algorithm framework, how could solution quality be improved from 95% to the exact optimum? For example, by extending annealing time, increasing iteration count, introducing multi-start parallel search, increasing the number of MTJs, or reducing system latency (to increase flip frequency)? If this is not achievable with the current framework, the authors should discuss possible strategies for progressing from near-optimal to exact solutions.\nAuthors’ response: We thank the reviewer for this important question. We believe that it remains an open challenge determining how to consistently progress from near-optimal to exact-optimal solutions in the QAP, and more generally in NP-hard optimization problems at scale. Nevertheless, several techniques can help to improve the solution quality.\nFrom an algorithmic perspective, the solution quality can often be improved by employing slower annealing schedules, longer annealing times (i.e., increased iteration counts), or more advanced annealing strategies such as parallel tempering, simulated quantum annealing, or adaptive temperature control. Algorithmic methods like Tabu search and Graver Augmented Multiseed Algorithm, will also enhance the search for feasible solutions. All these methods help the system escape local minima and explore the energy landscape more thoroughly. From a hardware perspective, increasing the number of MTJs allows larger-scale problem instances to be mapped directly onto the hardware, eliminating the overhead associated with MTJ reuse. Meanwhile, reducing system latency enables a higher number of sweeps (spin-flip attempts) within a given runtime. Both approaches can effectively enhance the overall solution quality.\nHowever, these heuristic improvements do not guarantee convergence to the exact global optimum. Due to the intrinsic computational complexity of the QAP, achieving exact solutions at scale remains an open and important research question. Nevertheless, in practical applications, the time-to-near-optimal solution and the energy required to reach it are often equally critical performance metrics. In this regard, probabilistic computing systems can play a pivotal role by accelerating the annealing process and significantly reducing the overall energy consumption required to approach near-optimal solutions efficiently. Such systems could be particularly well suited for complex optimization tasks in emerging edge and Internet of Things (IoT) applications, where energy efficiency and rapid decision-making are essential.\nIn response, we have added a discussion in the Main text (page 12) to provide readers with a more comprehensive understanding of possible approaches to achieving better solutions.\n(5) The manuscript reports a noticeable variation in MTJ switching probabilities (Fig. 2c) and compensates for it through calibration — a valuable approach to addressing device uniformity issues. From a device perspective, what might be the underlying causes of such variation?\nAuthors’ response: The observed variation in MTJ switching probabilities primarily arises from process-induced non-uniformities during fabrication. In particular, variations in thin-film thickness, and etching precision can lead to fluctuations in the device resistance, magnetic anisotropy, and MTJ pillar dimensions, all of which influence the switching probability. In response to the reviewer’s comment, we have added a brief explanation to the Main text (page 7).\nGiven the relatively narrow voltage tuning range and device-to-device inconsistency, could operational stability be affected (e.g., probability curve shifts over time or temperature, while stored σ and μ remain unchanged)?\nAuthors’ response: We thank the reviewer for this important point. Experimentally, we do observe a slight drift in switching probability with increasing device temperature, which can influence the stability of the probability curve if not properly managed. To mitigate this effect, we employ two practical measures:1. We precondition MTJs before operation to allow the system to reach thermal equilibrium, thereby minimizing transient thermal fluctuations; and 2. We use a small fan to prevent heat accumulation during prolonged measurements. These steps effectively stabilize the device temperature and ensure consistent switching behavior across experiments. In future implementations, an on-chip proportional-to-absolute-temperature (PTAT) compensation circuit could be integrated to automatically correct temperature-induced drift, thereby enhancing the long-term operational stability of the system.\nThe above discussion has been included in Supplementary Note 1 to provide readers with a more comprehensive understanding of the calibration process.\nCould changes in ambient magnetic fields impact global applicability?\nAuthors’ response: Ambient magnetic fields can influence MTJ characteristics, as the magnetization of the free and reference layers can be perturbed under strong external fields. However, commercial STT-MRAM chips have specifications regarding maximum external magnetic fields in the range of a few tens of mT to one hundred mT, which exceeds typical values found in everyday environments by two orders of magnitude [Dieny, B. et al., arXiv 2409.05584 (2024)]. This ensures that MTJ-based systems are globally applicable and reliable under normal operating conditions. For deployment in environments with substantial magnetic interference, further optimization of the MTJ stack [Xu, S. et al., Nat. Commun. 16, 8370 (2025)] and the implementation of appropriate magnetic shielding strategies [Dieny, B. et al., IEDM 2(3), 53–59 (2024)] may be necessary. Enhancing magnetic immunity and robustness has thus become an active area of research in spintronics.\nMoreover, the implications for large-scale scaling (e.g., millions of devices) are not discussed. Would storing σ and μ for all devices require significant FPGA memory, and would large-scale data transfer substantially increase computation time?\nAuthors’ response: We thank the reviewer for this insightful comment. As the core element of commercial magnetic random-access memory (MRAM) technology, the MTJ demonstrates outstanding device uniformity across large-scale integration. Therefore, we believe that for practical large-scale implementations, one set or several sets of σ and μ parameters can typically be sufficient to characterize the switching behavior of all MTJs, without consuming significant FPGA memory resources.\nWe agree with the reviewer that, as system scale increases, data transfer, especially between the memory and computing units can become a major performance bottleneck, thereby increasing the overall computation time. A direction to mitigate this issue could be to construct a hybrid computing architecture that integrates in-memory computing with probabilistic computing. For example, resistive random-access memory (RRAM) or magnetic random-access memory (MRAM)-based in-memory computing arrays could be combined with magnetic tunnel junction (MTJ)-based probabilistic computing units. Such a system would significantly reduce memory-computation data movement especially during MAC operations, enhance energy efficiency, and improve scalability for large-scale probabilistic Ising machines.\nIn response to the reviewer’s comment, we have added a discussion to the Main text (page 7), as “For large‐scale implementations (e.g., on the order of millions of MTJs), assigning distinct calibration parameters to every individual device would impose prohibitive hardware overhead. Importantly, wafer-level MRAM production has demonstrated excellent device-to-device uniformity in MTJ switching statistics. Based on this, it is reasonable to assume that only a small number of calibration parameter sets (e.g., one global set, or a few sets per process/wafer class) would be sufficient to capture the switching probability characteristics across all MTJs, without consuming significant memory resources. Furthermore, as the system scales, data transfer, particularly between memory and computing units, can substantially increase computation time. Benefiting from the CMOS compatibility, MTJ-based probabilistic computing systems can be integrated with emerging in-memory computing technologies such as MRAM51 and resistive random-access memory (RRAM)52 based platforms.”\n(6) In Supplementary Tables S1 and S2, the authors present power consumption values for FPGA, ADC+DAC, and MTJ components under different system scales. Were these values measured experimentally or calculated from device parameters and duty cycles? Providing more detail on the measurement or estimation methodology would enhance clarity.\nAuthors’ response: The power consumption of each system component was determined as follows. The FPGA power consumption was obtained from the Vivado Power Report generated after design implementation. The ADC and DAC power consumption values were calculated from their supply voltage and current specifications. The MTJ power consumption was estimated by adding the power consumption in the reset and perturb phase. P = I2(R +R ), where I (~50 µA) is the average switching current and R (~5k Ohm) and R (~2.5k Ohm) are the MTJ resistance in antiparallel state and parallel state, respectively. AP P AP P\nDetails on power consumption have been added to the Methods section and Supplementary Note 5.\n(7) On page 5, second paragraph, the full form of the abbreviation “PIM” is not given.\nAuthors’ response: Thank you for careful reading. We have added the full form of probabilistic Ising machine (PIM).\nResponse to the Second Reviewer\nReviewer #2 (Remarks to the Author):\nIn this work, the authors present an experimental demonstration of spin-transfer-torque magnetic tunnel junction (STT-MTJ) based probabilistic processor for solving quadratic assignment problems (QAPs). The integration of 144 MTJs with an FPGA represents a significant step forward in the field of spintronics-based probabilistic computing. The work is novel, well-executed, and the results are promising. However, several concerns must be addressed to strengthen the manuscript's clarity before it can be accepted for publication.\nAuthors’ response: We sincerely thank the reviewer for the supportive and constructive feedback, which has motivated us to improve our work. Please find our detailed, point-by-point responses below.\n(1) \"PIM\" (Probabilistic Ising Machine) is introduced without prior definition. Please introduce and define this term clearly for readers who may not be specialists in this specific sub-field.\nAuthors’ response: Thank you for careful reading. We have added the definition of probabilistic Ising machine (PIM).\n(2) Regarding Fig. 1d, the TSP energy landscape is depicted as smooth, while QAP is rugged. Is this true for all TSP instances, or just a selected case? The authors should clarify this point and provide appropriate citations to literature that discusses the comparative complexity of TSP vs. QAP landscapes to avoid over-generalization.\nAuthors’ response: We thank the reviewer for the comment. Fig. 1d was originally plotted from specific instances (a scale-4 QAP from TinyQAP dataset and TSP from the first four-city of Burma14 dataset) to showcase that QAP has more complex energy landscape compared to TSP. We agree with the reviewer that it is difficult to generalize the structural characteristics of high-dimensional combinatorial optimization landscapes from one specific small example.\nTo avoid over-generalization, we have replaced Fig. 1d with a schematic that instead highlights the localized interaction patterns underlying TSP, as shown in Fig. P4. The key intuition is that in the TSP, cities are connected in a chain where each city interacts only with its immediate neighbors along the route (local interactions), as shown in Fig. P5a, whereas in the QAP, every facility-location pair interacts with all others through both distance and flow terms, resulting in global all-to-all interactions, as shown in Fig. P5b. Specifically, TSP can be expressed as a special case of QAP. The TSP on n cities can be formulated as a QAP(F,D), where D is the city-city distance matrix, and F is the adjacency matrix of a Hamiltonian cycle on n vertices (F =1 if i and j are consecutive positions in the tour, 0 otherwise). The TSP cost then corresponds to the sum of distances between consecutive visiting cities. Thus, TSP corresponds to a QAP with a sparse F, whereas general QAP allows dense F. From a theoretical standpoint, TSP admits a polynomial-time 3/2-approximation [Christofides, N., Oper. Res. Forum 3, 20 (2022)], while QAP is not approximable in polynomial time within any finite factor unless P = NP [Queyranne, M., Oper. Res. Lett. 4, 231–234 (1986)]. The distinct interaction patterns lead to differences in computational complexity between the two problems. We believe that the revised figure and accompanying clarifications offer readers an intuitive yet rigorous understanding of why the QAP is inherently more challenging than the TSP.\nFig. P4| The TSP interaction pattern. Optimization considers only distances between consecutively visited cities.\nFig. P5| Structural differences between QAP and TSP. a, TSP structure, b, QAP structure.\n(3) Regarding the description of emerging spintronic technologies for probabilistic computing, it would be beneficial to include references on related works, such as: Nano Lett. 2024, 24, 5420−5428 (https://doi.org/10.1021/acs.nanolett.3c04820).\nAuthors’ response: We thank the reviewer for the suggestion. We have added the related reference.\n(4) The experiments are limited to 12-city QAPs, which may be too small to convincingly demonstrate scalability. What is the maximum scale (number of cities)? At what number of cities does this proposed system fail to achieve the optimal solution? Please test larger instances and report when the performance degrades or fails to find the optimum.\nAuthors’ response: We thank the reviewer for this comment. In the TinyQAP dataset employed for benchmarking, the largest problem instance has a scale of 12. QAP is known to be challenging, and the difficulty grows exponentially with problem scale. As shown in Fig. 5 of the main text, our system achieves optimal solutions consistently up to scale 7, while from scale 8 onward, the hardware begins to yield all near optimal but no optimal solutions, reflecting the intrinsic computational hardness of larger instances.\nMotivated by the reviewer’s advice, we performed additional scalability tests on larger Tai-series QAP instances from QAPLIB [Burkard, R. E., et.al, J. Global Optim. 10, 391–403 (1997)], and the results are summarized in Fig. P6. The simulated annealing (SA) scheme is implemented on MTJ-based probabilistic computing system, while for the parallel trial annealing (PA) scheme, it requires at least n2 parallel MTJs for a scale-of-n problem, which exceeds the capacity of our current proof-of-concept prototype. We perform simulations with the parallel-trial annealing scheme to show scalability. We tested both SA and PA under the same annealing schedule, where the inverse temperature increases linearly with respect to the sweep (for SA) or iteration (for PA). A sweep is defined as one full pass through all spins in the system, whereas an iteration corresponds to the update of a single spin. Thus, their time relation can be expressed as Number-of-sweep = N × Number-of-iteration, where N is the number of spins. Accordingly, since the SA process is annealed based on sweeps, its runtime is measured in sweeps, while the PA process is measured in iterations. In our experiments, we performed sweeps (equivalent to 2N×104 4 Markov-Chain-Monte-Carlo sampling iterations) for SA and iterations for PA. 2×10 5\nAs shown in Fig. P6a, for the simulated annealing (SA) scheme, the probability of finding feasible 2×10 solutions begins to decrease from 100% at a problem scale of 25, whereas for the parallel-trial annealing (PA) scheme, this decline starts at scale 40. In terms of solution quality, as shown in Fig. P6b, PA consistently achieves lower (better) normalized QAP energy compared to SA across all tested scales. It should be noted that under our hyperparameter settings, PA performs far fewer Markov Chain Monte Carlo (MCMC) sampling iterations than SA, particularly for large N. Nevertheless, PA outperforms SA in both feasibility probability and solution quality. With a larger number of iterations, PA is expected to achieve even better performance. Combined with the compact and energy efficient MTJ based p-bit implementation, MTJ-PA can be a strong candidate for scalable probabilistic computing schemes.\nThe above discussion has been added to Supplementary Note 2 to clarify the scaling characteristics of SA and PA.\nFig. P6 | Performance on Tai benchmark instances. a, Probability of feasibility as a function of QAP scale. b, The normalized QAP energy as a function of QAP scale. The normalized energy is defined as the ratio of the obtained QAP energy to the known optimum energy, where a value of 1 indicates an optimal solution.\n(5) Regarding Fig. 4c, the comparison is conducted between MTJ-PA and CPU-SA. Why was there no comparison with CPU-PA? To demonstrate the superiority of the MTJ-integrated system, a comparison between MTJ-PA and CPU-PA is essential.\nAuthors’ response: The key distinction lies in the parallel nature of the parallel-trial annealing (PA) algorithm. In a single-threaded CPU implementation, the parallel-trial annealing algorithm can only be emulated sequentially, as the CPU executes one instruction at a time. As a result, for a problem with n spins, the CPU based PA requires O(n) time per iteration to evaluate all spin flip proposals, which is effectively equivalent to one simulated annealing (SA) sweep. In contrast, on parallel hardware, all spins can be evaluated simultaneously, resulting in a time complexity of O(1) per iteration.\nAs shown in Fig. P7, we compare the time-to-solution of CPU-based SA, CPU-based PA, and MTJ-based SA and PA. For small problem sizes (e.g., scale below 7), the CPU-SA performs slightly better compared to CPU-PA because the CPU-PA implementation includes an additional arbiter process in each iteration to stochastically select among update candidates, which slightly increases the time consumption per iteration. However, starting from problem scale 8, CPU-PA begins to show modest improvement due to its ability to increase spin flip proposal acceptance probability. MTJ-based PA continues to outperform CPU-PA as the problem scales up, benefiting from the hardware-level parallelism inherent to the architecture, where all spin updates are evaluated concurrently.\nFig. P7| Time-to-solution comparison among CPU based PA, SA and MTJ based PA, SA.\n(6) The authors have mentioned “compared to CPU-SA, the MTJ-PA configuration achieves energy savings of up to 82%”. How was the energy consumption of CPU-SA measured? This manuscript does not mention how the energy/power consumption of CPU-SA was calculated/measured. Please provide a detailed methodology.\nAuthors’ response: The energy consumption of the CPU-SA was estimated based on measured power and runtime during execution. Specifically, the instantaneous CPU utilization of the MATLAB process was monitored using the Windows Task Manager while running the SA code. The effective CPU power was then obtained by multiplying the utilization percentage by the nominal power rating (125 W) of the CPU Intel i9-14900K. The execution time (time-to-solution) was recorded using MATLAB built-in tic and toc functions, which measure the wall-clock runtime of each experiment. The energy-to-solution was then computed as: E = P ×TTS, where P is the effective CPU power and TTS is the time-to-solution. We CPU CPU have added a detailed methodology to the Method section.\nResponse to the third Reviewer\nReviewer #3\nThis paper addresses an important topic and demonstrates interesting work, but requires revision. Following are my comments towards improving the manuscripts.\nAuthors’ response: We sincerely thank the reviewer for the constructive feedback. Please find our detailed, point-by-point responses below.\n1. Line 170: The estimate that an “MTJ-based implementation can reduce transistor count by over three orders of magnitude (~1200×) and energy consumption by more than 70× compared to conventional CMOS technology” is a powerful statement. However, it relies on a single citation (50). The validity of this comparison heavily depends on the assumptions made about the “conventional CMOS technology” in that reference. Briefly elaborate on the basis of this comparison. For instance, what kind of CMOS pseudo-random number generator (PRNG) and architecture is assumed, etc.\nAuthors’ response: We thank the reviewer for pointing out this. We agree that the term “conventional CMOS technology” was broad. The CMOS counterpart specifically refers to a CMOS-based pseudo-random number generator (PRNG) implemented using a 32-bit linear feedback shift register (LFSR). To realize the p-bit behavior (tunable probability), additional digital circuitry is required, specifically, look-up tables (LUTs) to store the hyperbolic tangent function and a 32-bit comparator to binarize the output, as illustrated in Fig. P8. We have revised our manuscript to include a more detailed description. Now the line reads: “It is estimated that an MTJ-based implementation can reduce transistor count by over three orders of magnitude (~1200×) and energy consumption by more than 70× compared to a conventional CMOS based 32-bit LFSR + LUT implementation53.”\nFig. P8| Schematic of CMOS based p-bit implementation.\n2. Figure 3: To improve clarity, consider adding a note in the caption, briefly explaining the arbiter’s function. This would help readers grasp its function without having to refer back to the main text.\nAuthors’ response: We thank the reviewer for the suggestion. We have improved the caption of Figure 3 by adding a brief introduction of the update arbiter block as “The update arbiter (highlighted in purple) uses the uniform random numbers generated by the MTJs to randomly select one spin from the set of candidate spins for update”.\n3. Thin film composition is provided in the “Methods” section, there is insufficient discussion of why this specific stack was chosen and how the perpendicular magnetic anisotropy was tuned for optimal stochastic switching.\nAuthors’ response: For large-scale probabilistic computing systems, device uniformity is of great importance. Significant device-to-device variation would introduce a heavy calibration burden, both in memory overhead (to store individual MTJ calibration parameters) and computational cost (to generate unique bias inputs). For this reason, we deliberately selected a conventional perpendicular STT-MTJ stack, which can offer better uniformity control, rather than superparamagnetic MTJs (with a thinner free layer thickness and smaller MTJ pillar size for smaller magnetic volume). This choice reflects a practical trade-off, higher energy consumption per switching in exchange for better device uniformity.\nAlthough our lab-fabricated MTJs exhibit some variability due to process-induced non-uniformity, as shown in Fig. P9 (details see Supplementary Note 1), we believe that commercial MRAM fabrication processes can achieve the uniformity required for scalable probabilistic computing systems, thereby positioning MRAM not only as a high-performance memory technology but also as a promising platform for probabilistic computing.\nFig. P9| Calibration parameter statistics of probabilistic switching responses of MTJs a, Distribution of slope factors. The slope factor determines the stretch of a sigmoid-like switching curve. b, Distribution of deviation factors. The deviation factor determines the horizontal deviation of the sigmoid-like switching curve from the standard sigmoid function which centered around zero.\n4. Figure 2c-d: Switching probability is mentioned but it lacks information on the non-linear device variations inherent in MTJ arrays. A more sophisticated calibration scheme accounting for thermal stability factor variations would be more appropriate.\nAuthors’ response: We thank the reviewer for the insightful comment. In brief, the slope factor defined in Supplementary Note 1 already captures the nonlinear device-to-device variation and is closely related to the thermal stability factor. Below, we provide a detailed analysis and show how the empirical slope factor in sigmoid-like function maps to the physics-based parameters (including thermal stability), thereby enabling a calibration scheme that accounts for these nonlinearities.\nAs discussed in Supplementary Note 1, we use an empirical sigmoid-like expression (1) 1 to fit each switching response, where p𝑃𝑃 i(s𝑉𝑉 t in h)e =switchin𝜎𝜎g( 𝑉𝑉pirno−b𝜇𝜇a)bility, σ is the slope factor (non-linear), µ is 1+𝑒𝑒 the deviation factor (linear), and V is the input voltage to the p-bit cell. From device physics, STT-MTJ in switching probability can be expressed as [Bedau, D. et al., Appl. Phys. Lett. 97, 262502 (2010)]: , (2) 𝑡𝑡 𝐼𝐼 where t is the duration of the in𝑃𝑃pu=t c1ur−renetx ppul−se,𝜏𝜏 0e xips th−e Δatte1m−pt𝐼𝐼 𝑐𝑐t0ime, is the thermal-stability parameter, and is the critical switching current. In this s 𝜏𝜏 tu0dy, the MTJs are bia Δ sed with voltage pulses, therefore, we s 𝐼𝐼 ubstitute I and I c0 with V and V c0 respectively, without changing the mathematical structure: . (3) 𝑡𝑡 𝑉𝑉 Linearizing eq. (3) at the midp𝑃𝑃oin=t P1 =− 0e.x5p giv−e𝜏𝜏s0 ae lxopca−l sΔlop1e:− 𝑉𝑉𝑐𝑐0 . (4) 𝑑𝑑𝑑𝑑 Δ𝑙𝑙𝑙𝑙2 The midpoint slop of the sigmoid-like |fu𝑑𝑑n𝑉𝑉c|t𝑑𝑑io=n0 .5is= 2𝑉𝑉𝑐𝑐0 . For a specific MTJ probabilistic switching response, we can equate these two slope values, which gives −𝜎𝜎/4 , (5) −2Δ𝑙𝑙𝑙𝑙2 indicating the relation between the slope fact𝜎𝜎or= an𝑉𝑉d𝑐𝑐0 thermal stability.\nFig. P10a and Fig. P10b. show that both the empirical model (Fig. P10a) and the physics-based 𝜎𝜎 model (Fig. P10b) fit the switching data with high accuracy. Fig. P10c and P10d present the distributions of the sigmoid function slope factor and the physics-based thermal stability factor Δ. As derived in Eq. (5), these parameters are quantitatively related.\nAlthough both methods accurately fit the experimental data, the empirical model Eq. (1) provides a more intuitive link between the algorithm and the physical input voltage amplitude. The input to p-bit I 19 can be simply expressed as: . This expression makes it more directly applicable for system- 𝜎𝜎(𝑉𝑉in−𝜇𝜇) level implementation. Thus, we adopt the empirical method for our system setup. 𝐼𝐼 = 𝑒𝑒\nFig. P10| Comparison of empirical and physics-based calibration methods. MTJ probabilistic switching responses data and fitting lines from empirical method (a), and physics-based method (b). Distribution of slope factor (c), and thermal stability factor (d).\n5. Figure 4: The energy consumption analysis includes FPGA, DAC, and ADC power but lacks breakdown of individual components. This makes it difficult to assess where optimization efforts should focus. The 93% energy savings claim requires more detailed power budgeting.\nAuthors’ response: We appreciate this comment. A detailed power breakdown of each system component has been provided in Tables S1 and S2, including contributions from the FPGA, DACs, and ADCs. Among these components, the FPGA accounts for the largest portion of power consumption (> 90% in MTJ-SA and > 70% in MTJ-PA), primarily due to DAC/ADC control and multiply-accumulate (MAC) operations. The off-chip commercial DACs and ADCs also contribute ~ 5% in MTJ-SA, and ~30% in MTJ-PA, as they are not optimized for low-power operation. To further enhance energy efficiency, the breakdown in Tables S1 and S2 indicates that future optimization efforts could focus on integrating low-power DACs, ADCs, and control logic on-chip, which would substantially reduce FPGA control and off-chip interface power, thereby improving overall system efficiency.\nResponse to the Fourth Reviewer\nReviewer #4 (Remarks to the Author):\nClassical probabilistic computing, in particular for topics of monte carlo sampling, is certainly of interest in terms of computational technologies. The comparison to simulated annealing is good to show. The comparison against D-Wave quantum annealers is a somewhat strange comparison choice in my view -- typically, especially for these types of optimization problems, simulated annealing performs much better than D-Wave quantum annealers. Most importantly, having the ground state optimal solution found by Gurobi is an excellent part of the study -- this is crucial to show if the evaluated methods can reach the ground-state energy.\nAuthors’ response: We thank the reviewer for recognizing the relevance of classical probabilistic computing and for appreciating our comparison with simulated annealing. Our comparison with D-Wave quantum annealers is to offer a hardware benchmark across digital, probabilistic, and quantum annealing platforms. We agree that, for many combinatorial optimization problems, classical simulated annealing often achieves better performance than current D-Wave systems due to embedding overhead and device noise. The D-Wave results are included to contextualize our MTJ-based probabilistic annealer within the broader family of annealing hardware (digital, probabilistic and quantum). We also greatly appreciate the positive remark on the Gurobi baseline.\nSpecific comments:\n1. The methods used in the D-Wave QPU runs needs to be specified in complete detail. Namely, what minor embeddings were used, chain strength, broken chain repair algorithm. Also, importantly the total number of anneals measured needs to be stated (and whether the solution quality is an average, or the best). Also, the total number of logical variables in the original QUBO needs to be stated for complete methods.\nAuthors’ response: We thank the reviewer for this comment. The details of the D-Wave experiments, including the minor-embedding method, chain strength, broken-chain repair algorithm, number of anneals, and number of logical variables in the original QUBO, have been summarized in Table P2 (and also included in Supplementary Note 6). To evaluate the solution quality, 100 anneals are performed for each instance, and all feasible solutions are collected to calculate the mean and standard deviation, as shown in Fig. 5 of the main text and described in its caption.\nImplementationa Advantage_system6.4 Advantage2_system1.6\nGraph topology Pegasus Zephyr\nMinor embedding Minorminer (default) Minorminer (default)\nmethod\nChain strengthb Scaled to the problem (default) Scaled to the problem (default)\nBroken chain repair Majority vote (default) Majority vote (default)\nalgorithm\nAnnealing time 2000 µs 2000 µs\nNumber of anneals 100 100\nQUBO problem 9,16,25,36,49,64,81,100,121, and 9,16,25,36,49,64,81,100,121, and\nscale 144 144\naAs Advantage_system7.1 and Advantage2_prototype2.6 systems are not currently available, all testing results are updated with the latest system Advantage_system6.4 and Advantage2_system1.6.\nbDetailed chain strength of each testing instance is listed in Table. P3 and P4.\n2. Figure 1d: it should be clarified that this is purely a schematic, not actual data, it is somewhat confusing otherwise. Line 93 on page 5 states that \"By contrast, the TSP plot looks much more like a smooth bowl\", which is confusing to me -- is this an empirical statement? I think it is somewhat hard to correctly interpret the search space of a high dimensional discrete combinatorial optimization search space. This portion of the manuscript needs to be clarified as to whether this is actual data, or is a conceptual interpretation of the search spaces.\nAuthors’ response: We thank the reviewer for raising the point. Figure 1d was originally generated from specific instances (a scale-4 QAP from TinyQAP and TSP from the first four-city of Burma14) to showcase that all-to-all connected QAP has more complex energy landscape compared to TSP. We agree with the reviewer that it is difficult to correctly interpret the structural characteristics of high-dimensional combinatorial optimization landscapes.\nTo avoid over-generalization, we have replaced Fig. 1d with a schematic that instead highlights the localized interaction patterns underlying TSP, as shown in Fig. P11. The key intuition is that in the TSP, cities are connected in a chain where each city interacts only with its immediate neighbors along the route (local interactions), as shown in Fig. P12a, whereas in the QAP, every facility-location pair interacts with all others through both distance and flow terms, resulting in global all-to-all interactions, as shown in Fig. P12b. Specifically, TSP can be expressed as a special case of QAP. The TSP on n cities can be formulated as a QAP(F,D), where D is the city-city distance matrix, and F is the adjacency matrix of a Hamiltonian cycle on n vertices (F =1 if i and j are consecutive positions in the tour, 0 otherwise). The TSP cost then corresponds to the sum of distances between consecutive visiting cities. Thus, TSP corresponds to a QAP with a sparse F, whereas general QAP allows dense F. From a theoretical standpoint, TSP admits a polynomial-time 3/2-approximation [Christofides, N., Oper. Res. Forum 3, 20 (2022)], while QAP is not approximable in polynomial time within any finite factor unless P = NP [Queyranne, M., Oper. Res. Lett. 4, 231–234 (1986)]. The distinct interaction patterns lead to differences in computational complexity between the two problems. We believe that the revised figure and accompanying clarifications offer readers an intuitive yet rigorous understanding of why the QAP is inherently more challenging than the TSP.\nFig. P11| The TSP interaction pattern. Optimization considers only distances between consecutively visited cities.\nFig. P12| Structural differences between QAP and TSP. a, TSP structure, b, QAP structure.\n3. I would generally say that the methods description, and properties, of the D-Wave quantum annealing devices is not very sufficient -- it would be good to cite and describe the technology ideas, how it works, why it may be reasonable to use in this context, etc. Maybe also why minor embedding is required, maybe also how much of the QPU chip was used in the embedding process.\nAuthors’ response: We thank the reviewer for the suggestion. We have added a paragraph in the Main text describing the operating principles of the D-Wave quantum annealer and the justification for minor embedding as follows.\n“The D-Wave annealer employs a process known as quantum annealing to solve optimization problems by finding the lowest-energy configuration of a system of superconducting qubits. Each optimization problem is first formulated as an energy landscape, defined by setting the qubit biases and coupler strengths that represent the objective function of the problem. The system begins with all qubits in a quantum superposition of possible states. During the annealing process, the Hamiltonian parameters are gradually evolved so that quantum fluctuations and entanglement guide the system through many possible configurations. As the annealing progresses, the influence of the driver Hamiltonian is reduced while that of the problem Hamiltonian is increased, allowing the system to evolve naturally toward low-energy states [Johnson, M. W. et al., Nature 473, 194–198 (2011)]. Because the connectivity of the D-Wave hardware graph (e.g., Pegasus and Zephyr topologies) does not always match the connectivity of the problem graph, minor embedding is required to map logical variables and couplings onto the available physical qubits. In this process, a single logical variable may be represented by a chain of physically connected qubits, with the chain strength tuned to preserve variable consistency during annealing. At the end of the annealing, typically within time scale of microseconds, the quantum superposition collapses, and the qubits settle into a single classical configuration corresponding to the minimum (or near-minimum) energy of the encoded problem, representing the final solution. As a commercially available quantum computing platform, D-Wave has demonstrated the ability to address certain real-world optimization problems. In this work, we use D-Wave as a representative quantum annealer to provide a hardware benchmark alongside the digital CPU and our MTJ-based probabilistic annealer.”\nDetailed QPU chip utilization rates for both hardware topologies across all QAP instances are summarized in Tables P3 (Pegasus) and P4 (Zephyr). The number of utilized physical qubits can be estimated by multiplying the physical qubit usage percentage by the total number of qubits. These two tables are added to Supplementary Note 6.\nInstances Logical Logical Physical Physical Avg. chain Chain\nvars. couplers qubits usage couplers length strength\nusage\nTiny03a 9 36 0.25% 0.12% 1.56 1.25\nTiny04a 16 120 0.71% 0.44% 2.5 1.99\nTiny05a 25 300 1.53% 1.09% 3.44 3.54\nTiny06a 36 630 3.26% 2.46% 5.08 4.40\nTiny07a 49 1176 6.22% 5.31% 7.12 6.19\nTiny08a 64 2016 9.41% 8.44% 8.25 8.07\nTiny09a 81 3240 15.86% 14.24% 10.99 11.04\nTiny10a 100 4950 22.08% 20.21% 12.39 14.24\nTiny11a 121 7260 32.34% 29.67% 15.00 16.19\nTiny12a 144 10296 54.90% 49.36% 21.40 17.13\nTable. P4| QPU (Zephyr) chip utilization rates. Total number of qubits is 4593, and total number of couplers 41769.\nInstances Logical Logical Physical Physical Avg. chain Chain\nvars. couplers qubits usage couplers length strength\nusage\nTiny03a 9 36 0.3% 0.13% 1.56 1.25\nTiny04a 16 120 0.72% 0.42% 2.06 1.99\nTiny05a 25 300 1.46% 1.01% 2.68 3.54\nTiny06a 36 630 3.24% 2.40% 4.14 4.40\nTiny07a 49 1176 6.73% 5.17% 6.31 6.19\nTiny08a 64 2016 9.36% 7.95% 6.72 8.07\nTiny09a 81 3240 16.37% 14.43% 9.28 11.04\nTiny10a 100 4950 26.52% 24.30% 12.18 14.24\nTiny11a 121 7260 49.31% 43.95% 18.72 16.19\nTiny12a 144 10296 58.87% 53.62% 18.78 17.12\n25\n4. The CPU compute time used by Gurobi should be reported as well.\nAuthors’ response: The below Gurobi implementation information has been included to Supplementary Note 8.\nCPU model: INTEL(R) XEON(R) GOLD 6526Y\nGurobi 12.0.3\nTable. P5| Gurobi runtime for solving QAP instances.\nInstance Problem scale QAP energy Time-to-traverse (s) Time-to-\nsolution (s)\nTiny03a 3 2.34 0.0314 0.0216\nTiny04a 4 2.21 0.0383 0.0274\nTiny05a 5 5.72 0.1082 0.0385\nTiny06a 6 6.78 0.1748 0.1357\nTiny07a 7 9.35 0.5083 0.3167\nTiny08a 8 9.83 0.9375 0.6652\nTiny09a 9 16.26 2.6897 2.0785\nTiny10a 10 20.27 19.9439 5.2512\nTiny11a 11 22.72 104.6843 73.3535\nTiny12a 12 29.23 596.5185 147.7377\n5. One overall comment I have on the aims and claims of this study is that in order to really make a state of the art methods comparison between a hardware probabilistic sampler and classical digital computers would be to also run a GPU comparison -- highly optimized QPU simulated annealing implementations can be much more efficient than CPU implementations.\nAuthors’ response: We thank the reviewer for raising this important point. In general, GPUs can accelerate simulated annealing (SA) mainly through two mechanisms: (i) accelerating ΔE evaluations, and (ii) running a large number of SA trajectories in parallel and selecting the best solution among them. Motivated by the reviewer’s suggestion, we implement a GPU-optimized SA on an NVIDIA H100 (C++/CUDA) and benchmark it on the same TinyQAP dataset.\nBefore benchmarking GPU performance with the other computing schemes (CPU-SA, MTJ-SA, and MTJ-PA), we first focused on optimizing the GPU implementation itself, specifically along the two acceleration mechanisms discussed above.\n26\n1. ΔE evaluation: To ensure a fair comparison, the GPU implementation also uses the same incremental ΔE update strategy as the CPU baseline, rather than recomputing two full matrix-vector products at every iteration. (some details: As there is at most one spin update in each iteration in our simulated annealing configuration, the change in energy ( where H is the Ising Hamiltonian) can be expressed as , where 𝛥𝛥𝐸𝐸𝑖𝑖 = 𝐻𝐻|−𝑠𝑠𝑖𝑖, −wi𝐻𝐻th| 𝑠𝑠s𝑖𝑖 ,denoting the spin state, J the coupling matrix and h the loc 𝛥𝛥 𝐸𝐸 = l 𝑖𝑖fie ld − te 2 r 𝑠𝑠 𝑚𝑚𝑖𝑖𝐿𝐿 . 𝐿𝐿 𝑖𝑖If the ini 𝐿𝐿 ti𝑖𝑖a 𝐿𝐿 l L∑ i𝑖𝑖s 𝐽𝐽 p𝑖𝑖𝑖𝑖r 𝑠𝑠 e𝑖𝑖-c al ℎ l𝑖𝑖ated, then in each iteration, L can be updated incrementally as , where , if previous spin (s ) update i-1 proposal is accepted, and 𝐿𝐿 ot𝑖𝑖h = erw 𝐿𝐿𝑖𝑖is + e. 𝛿𝛿 I 𝐿𝐿 n𝑖𝑖 this case 𝛿𝛿 , th𝑖𝑖e = H − a 2 m 𝑠𝑠 i𝑖𝑖l−to1n 𝐽𝐽𝑖𝑖i,a𝑖𝑖−n1 update is effectively reduced to one addition operation, tw 𝛿𝛿 o 𝐿𝐿 c𝑖𝑖o = nd 0 itional (if-else) evaluations, and two scalar multiplications.)\nGPU acceleration is achieved by parallelizing the computation of the vector L at each iteration. In contrast to the CPU implementation, where the elements of L are updated sequentially, the GPU computes all n elements in parallel. Consequently, the effective per-iteration complexity for updating L is reduced from O(n) on the CPU to O(1) on the GPU (assuming sufficient parallel resources for an n-spin problem).\n2. Parallel replicas: Running many independent SA trajectories in parallel increases the probability of discovering low-energy solutions within a given wall-clock time. Since each trajectory explores a different stochastic search path, selecting the best solution among a large ensemble of replicas improves both expected solution quality and time-to-solution.\nSpecifically, we implement the ΔE update computations within CUDA threads, and we map parallel SA replicas onto separate CUDA thread blocks. We then conduct two systematic studies: (i) to determine the optimal number of replicas, and (ii) to verify that the H100 hardware resources can support fully parallel ΔE evaluations for each test instance.\nIn Fig. P13a, we examine how the time-per-sweep scales with the number of replicas executed on the GPU. Here, “time-per-sweep” denotes the wall-clock time required to complete one simulated-annealing sweep (i.e., to perform n spin-flip proposals, where n is the total number of spins). We observe that the time-per-sweep begins to increase sharply beyond 1024 replicas, indicating that the available GPU parallel resources are saturated, and additional replicas are therefore forced to execute partially sequentially due to hardware scheduling limits. In Fig. P13b, we analyze the scaling of time-per-iteration, where “time-per-iteration” refers to the wall-clock time required to complete a single spin-flip proposal evaluation. We observe that, for all QAP instance sizes (number of spins ranging from 9 to 144), the time-per-iteration remains essentially constant when the number of replicas is below 1024. This indicates that the per-iteration cost does not grow with the number of spins, i.e., the effective complexity is O(1) rather than O(n). Consequently, we confirm that the H100 hardware resources are sufficient to support fully parallel ΔE evaluations for all test instances when operating with fewer than 1024 replicas.\nWe next investigate the algorithmic benefit obtained from running replicated SA trajectories, as summarized in Fig. P13c. Here, we use sweep-to-solution as the figure of merit, defined as the number of 95 sweeps required to reach 95% of the optimal solution. We observe that increasing the number of replicas reduces the sweep-to-solution, indicating that parallel replication improves the probability of rapidly reaching high-quality solutions. However, trade-offs should be made between the sweep-to-solution benefits and the time-per-sweep overhead caused by limited hardware resources above 1024 replicas. Fig. P13d presents the scaling of TTS95 with respect to the number of replicas. Beyond 1024 replicas, especially for larger instances (QAP scale > 7), the algorithmic benefits are not sufficient to compensate for the increased time-per-sweep overhead. Consequently, we choose 1024 replicas as the optimal configuration to obtain the most efficient SA performance.\nFig. P13| GPU based replicated simulated annealing optimization. Time per sweep (a), time per iteration (b), sweep to solution (c) and TTS (d) scaling with regard to number of replicas for different 95 95 instance scales.\nFig. P14 compares the MTJ-based probabilistic sampler with digital samplers, including CPU and GPU. The GPU implementation uses the optimized configuration discussed above (1,024 replicas with fully 28 parallel L computation in each replica). The CPU implementation is based on an Intel i9-14900K processor, evaluated using both a compiled C++ implementation and an interpreted MATLAB implementation. Fig. P14a shows the comparison in terms of time-to-solution, where GPU performance is broadly comparable to CPU performance. Fig. P14b presents the energy-to-solution comparison. For the GPU-based replicated SA, we estimate a baseline power consumption of 90 W for the NVIDIA H100. We observe that the GPU consumes approximately two orders of magnitude more energy than the CPU.\nFig. P14| Performance comparison between MTJ based probabilistic sampler and digital samplers, including CPU and GPU. a, Time to solution comparison. b, Energy to solution comparison.\nThe fact that GPU-replicated SA achieves similar time-to-solution as the CPU while consuming more energy under the current benchmark instances is reasonable. As shown in Fig. P15, the time-per-iteration on the GPU is roughly 4 to 5× slower than on the CPU, even though the GPU computes L in parallel. This means that, in order for the GPU-replicated SA to ultimately outperform the CPU in time-to-solution, it would need to reach the target solution in only ~1/4 to 1/5 as many sweeps. However, as shown in Fig. P13c, particularly for larger instances, this reduction in sweeps is not achieved. To demonstrate a GPU advantage, the problem size must be large enough that the cost of evaluating the local energy change ΔE (computing the local field L) dominates over parallelization overheads, otherwise, CPU implementation can match or even outperform the GPU. This behavior is consistent with GPU performance studies, which report that CPUs often perform better for small problem sizes, whereas GPUs only show strong speedups once the workload is sufficiently large and highly parallel [Lee, V. W. et al., Proc. Int. Symp. Comput. Archit., 451–460 (2010)].\nFollowing the reviewer’s suggestion, we now include the GPU-based replicated SA performance comparison in the Main Text. The associated optimization details are provided in Supplementary Note 3.\n6. Line 243: \"For digital computing, we implement a simulated annealing algorithm on a CPU with 5000 sweeps\" For this type of computing paradigm comparison, the exact implementation really matters. Therefore, in this case, I think the authors should give more details on the programming language used, ideally a specific library that was used. Ideally this would be a compiled language, something along the lines of C, Julia, Rust, in order for the comparison to be very rigorous.\nAuthors’ response: We thank the reviewer for the suggestion. Our CPU implementation is based on MATLAB programming (line 209), which is commonly used in p-bit research field [Aadit, N.A. et al., Nat. Electron. 5, 460–468 (2022)], [Singh, N.S. et al., Nat. Commun. 15, 2685 (2024)], [He, Y. et al., IEEE J. Explor. Solid-State Comput. Devices Circuits 9, 83–91 (2023)]. While the choice of programming language does not affect the algorithmic behavior, the use of an interpreted language increases the runtime compared to a compiled one, due to additional overhead in code execution. We agree with the reviewer that implementing the algorithm in a compiled language would provide a more rigorous benchmark of the absolute CPU performance.\nMotivated by the reviewer’s suggestion, we implemented the same simulated annealing algorithm in C++. To ensure a fair comparison and avoid the influence of matrix multiplication performance differences between MATLAB and C++, we adopted an algorithmic form that does not involve direct matrix-vector multiplication between the J matrix and the spin vector. As there is at most one spin update in each iteration in our simulated annealing configuration, the change in energy ( where H is the Ising Hamiltonian) can be expressed as , where 𝛥𝛥𝐸𝐸𝑖𝑖 = 𝐻𝐻|−, 𝑠𝑠w𝑖𝑖 −ith𝐻𝐻 s| 𝑠𝑠d𝑖𝑖e,noting the spin state, J the coupling matrix and h the loca 𝛥𝛥 l 𝐸𝐸 f𝑖𝑖ie = ld − te 2 rm 𝑠𝑠𝑖𝑖. 𝐿𝐿 𝑖𝑖If the init 𝐿𝐿 ia𝑖𝑖l = L i∑s 𝑖𝑖p 𝐽𝐽 p𝑖𝑖e𝑖𝑖- 𝑠𝑠 c𝑖𝑖a lcu ℎ l𝑖𝑖ated, then in each iteration, L can be updated incrementally as , where , if previous spin (s i i- 1 ) update proposal is accepted, and 𝐿𝐿 ot𝑖𝑖h = erw 𝐿𝐿𝑖𝑖is + e. 𝛿𝛿 I 𝐿𝐿 n𝑖𝑖 this case 𝛿𝛿 , th𝑖𝑖e = H − a 2 m 𝑠𝑠 i𝑖𝑖l−to1n 𝐽𝐽𝑖𝑖i,a𝑖𝑖−n1 update is effectively 𝛿𝛿𝐿𝐿𝑖𝑖 = 0\n30\nreduced to one addition operation, two conditional (if-else) evaluations, and two scalar multiplications, thereby bypassing the influence of matrix–vector multiplication performance differences between MATLAB and C++.\nFurthermore, we have improved our MTJ-based probabilistic computing system. In our previous MTJ based parallel-trial (PA) scheme, a key speed bottleneck originated from the number of MTJs. For a QAP instance of size n, we previously employed n² parallel MTJ cells. However, this configuration required reusing MTJs within each iteration to generate random numbers for the Markov-Chain-Monte-Carlo (MCMC) sampling block and the update arbiter block, which introduced significant latency. We now integrate more MTJs (2n² – 1 for an instance of scale n). This eliminates the need to reuse MTJs during each iteration, thereby unleashing the full potential of our parallel MTJ-based probabilistic computing architecture. To realize this, we used additional MTJs and a second MTJ carrier printed circuit board. The two MTJ boards are connected to the field-programmable gate array (FPGA) via IPEX cables and can be controlled in parallel by the FPGA.\nThe benchmarking results among MTJ-based PA, MTJ-based SA, and CPU-based SA implemented in interpreted (MATLAB) and compiled (C++) languages for QAP instances of increasing scale are shown in Fig. P16. As seen from Fig. P16a, under identical annealing schedule, CPU-SA implemented in C++ outperforms that in MATLAB, confirming that compiled language implementations provide higher computational efficiency than interpreted ones. Furthermore, the corresponding energy-to-solution analysis is presented in Fig. P16b. The energy-to-solution is obtained by multiplying the system power by the corresponding time-to-solution. For the CPU-based implementation, the power consumption is estimated by multiplying the CPU utilization rate of the program (MATLAB or C++ executable) during runtime by the nominal power rating of the Intel i9-14900K processor. For the MTJ based implementation, the total power consumption is estimated by summing the power of FPGA (from Vivado implementation power report), DACs, ADCs (specified in their datasheets), and MTJs (calculated from input voltage and device resistance). The acceleration achieved by MTJ-PA compared to other implementations is defined as TTS TTS Compared with MTJ-SA, MTJ-PA achieves an acceleration factor of up to 123. other scheme MTJ-PA\nCompared wit / h CPU-base.d SA, MTJ-PA achieves up to 6.6× acceleration for the interpreted-language implementation and 3.2× for the compiled-language implementation. The energy saving factor is defined as ETS ETS . MTJ-PA achieves up to 98.3% energy reduction compared to MTJ- MTJ-PA other scheme\nSA, 87.1% compared to the interpreted-language CPU-SA, and 58.3% compared to the compiled-language 1−( / ) CPU-SA. It should be noted that our MTJ based probabilistic computing system is a proof-of-concept prototype with considerable potential for optimization, such as the integration of higher-speed/low power DACs and sense amplifiers for MTJ biasing and readout. We have incorporated the above discussion into the Main text.\n7. Line 173: in regards to the random number generation, this point is interesting because thermodynamic based random number generation is interesting. However, it is not clear in this study what exactly is the advantage gained by this particular random number generator. Generally, I would argue that a dedicated analysis of the random number generator would be good (e.g. min entropy, and general random number generator tests such as from Diehard and the NIST testsuite) -- and then, showing what exactly is \"high-quality\" (to quote the manuscript) about this randomness.\nAuthors’ response: We thank the reviewer for the comment. The efficient realization of a tunable random number generator is a key advantage gained by the MTJ-based implementation. In this work, STT-MTJs are employed not merely as random number generators for randomly selecting spins to update, but as tunable random number generators (p-bits in the context of probabilistic computing), where the probability of generating 0s or 1s can be flexibly adjusted by the input bias. This tunability is the essential advantage of MTJ-based random number generators in our system.\nWhen comparing the CMOS based p-bits and MTJ based p-bits. CMOS p-bit implementations typically rely on a pseudo-random number generator (e.g., linear feedback shift registers, LFSRs) combined with a look-up tables (LUTs) to achieve non-linear tunable probability distribution, which incurs additional hardware overhead. The MTJ-based p-bit, by directly exploiting thermally activated stochastic switching, 32 provides a compact, inherently probabilistic, and scalable solution. As reported by Singh et al. [Nat. Commun. 15, 2685 (2024)], MTJ-based implementations can reduce transistor count by more than three orders of magnitude (~1200×) and energy consumption by over 70× compared with CMOS counterparts (LFSR+LUT). Consequently, the key advantage of the MTJ-based random number generator lies in its ability to provide parallel, high-speed, and tunable random sources, which are essential for probabilistic sampling processes that random numbers are extensively required.\nIn the Main text, our intention in using the term ‘high-quality’ randomness was to contrast the MTJ-based random number generator with the linear feedback shift register (LFSR)-based pseudo-random generator, whose output exhibits linear dependency and periodic repetition after a certain sequence length, characteristics we consider of “low-quality” randomness. In contrast, the MTJ-based random number generator produces stochastic outputs that do not exhibit linear dependency or periodicity, regardless of runtime, owing to its thermally activated switching nature. However, we agree with the reviewer that the term ‘high-quality’ can be vague in this context, as from a statistical perspective, similar to other true random number generator implementations, the raw output of our MTJ-based source does not fully satisfy stringent statistical tests such as the NIST 800-22 suite without post-processing. As demonstrated in our prior work ([Ng et al., Phys. Rev. Appl. 19, 034077 (2023)], the NIST 800-22 results are summarized in Table P4 under “Raw data”), the unprocessed sequence does not pass all the tests. Furthermore, analysis using the NIST 800-90B indicates a min-entropy of 0.778584 bits per output bit, confirming that the entropy is below the ideal value of 1.\nFor security-oriented applications, such entropy reduction necessitates post-processing to enhance statistical uniformity. Accordingly, we applied Toeplitz hashing to the raw output, and the processed (Extracted) data successfully passed the NIST SP 800-22 tests, as shown in Table P6. Thus, in this case, the MTJ should be regarded as an entropy source for a true random number generator, which requires an additional post-processing unit to form a complete TRNG system.\nHowever, for probabilistic computing tasks such as MCMC sampling, the algorithmic performance is not highly sensitive to exact statistical distributions. Instead, scalability, tunability, and speed are the more dominant factors. Therefore, in this study, post-processing for statistical enhancement was not applied, as the focus was on the tunable probabilistic behavior.\nIn the Main text, we have therefore replaced the term “high-quality randomness” with “inherent and tunable randomness.” The revised sentence now reads: “The inherent and tunable randomness enhances the efficiency of the system to escape local minima and explore the energy landscape more effectively.”.\nRaw data Extracted data\nStatistical test\nP-value Proportion Result P-value Proportion Result\nFrequency 0.0000 0.1901 Failure 0.2315 1.0000 Success\nBlock frequency 0.0000 0.8651 Failure 0.9233 0.9882 Success\nCumulative sums 0.0000 0.1839 Failure 0.4201 1.0000 Success\nRuns \\ \\ \\ 0.8606 0.9941 Success\nLongest run 0.0000 0.1797 Failure 0.6436 1.0000 Success\nRank 0.1806 1.0000 Success 0.4547 0.9822 Success\nFFT 0.0000 0.5253 Failure 0.4666 1.0000 Success\nNonoverlapping\n0.0000 0.0000 Failure 0.0001 0.9704 Success\ntemplate\nOverlapping\n0.0000 0.0335 Failure 0.0024 0.9941 Success\ntemplate\nUniversal 0.0000 0.0000 Failure 0.6825 0.9882 Success\nApproximate\n\\ \\ \\ 0.5659 0.9941 Success\nentropy\nRandom excursions 0.0003 0.8182 Failure 0.2307 0.9703 Success\nRandom excursions\n0.0909 0.9091 Success 0.0780 0.9604 Success\nvariant\nSerial 0.0000 0.0000 Failure 0.1597 0.9763 Success\nLinear complexity 0.3781 0.9862 Success 0.3979 0.9882 Success\n8. Selection of the annealing time on the D-Wave QPUs: a choice if 2000 microseconds is reasonable, generally it should result in the lowest sampled energy. However, it is possible for some types of optimization problems, a slightly shorter annealing time is better. Ideally, one would perform some sort of systematic annealing time analysis. But, I do not think this crucial for the article to be published.\nAuthors’ response: Motivated by the reviewer’s comment, we systematically examined the influence of annealing time on each QAP instance. Specifically, annealing times of 2, 20, 200, and 2000 µs were tested on both Pegasus and Zephyr topologies, with 100 anneals per condition to obtain mean and standard deviation values. As shown in Fig. P5, the normalized QAP energy and feasibility probability were evaluated as a function of problem scale. As the reviewer pointed out, for certain instances, shorter annealing times can yield better results. For example, in the Pegasus topology (Fig. P17a), 20 µs outperforms 200 and 2000 µs across scales 3 to 5. In contrast, the Zephyr topology (Fig. P17c) shows that 200 µs performs best at scales 3 and 4, while 2000 µs performs best at scale 5.\nBoth Pegasus and Zephyr topologies fail to find feasible solutions beyond scale 5; thus, data points in Fig. P17a,c terminate at this scale. To assess overall performance, we further analyzed the Ising energy (the energy from the Ising Hamiltonian) without imposing feasibility constraints, as shown in Fig. P18.\nFig. P17| Influence of annealing time on D-Wave performance across problem scales and topologies. (a,b) Results obtained using the Pegasus topology and (c,d) using the Zephyr topology. (a,c) Normalized QAP energy and (b,d) probability of feasibility as a function of QAP scale for annealing times of 2, 20, 200, and 2000 µs. Each data point represents the mean value from 100 anneals, with error bars indicating standard deviations. The dashed line denotes the optimal energy.\n9. Line 80, page 4: \"function becomes equivalent to finding the ground state spin configuration\": is there any potential degeneracy in this formulation? This sentence implies that there is only ever one unique solution.\nOverall I would not recommend this manuscript for publication in its current form -- it requires revisions to clarify the methods, and some of the plots. In terms of general community interest, I would estimate that this study is not of considerable interest to the quantum computing community. But, I do think the most interesting part of the study is the hardware description of the probabilistic Ising model sampler.\nAuthors’ response: As the reviewer pointed out, the main focus of our work is the experimental demonstration of the MTJ based probabilistic Ising model sampler. Our discussion on digital and quantum computing is intended as a benchmarking reference to contextualize our MTJ-based probabilistic annealer in the broader family of annealing hardware, including digital, probabilistic and quantum. We believe such cross-paradigm comparisons can help readers gain a more comprehensive understanding of the probabilistic computing scheme. We have revised the manuscript to clarify the methods and improve the plots and analysis according to reviewer’s constructive comments. We believe that the changes resolve the earlier concerns, as well as presenting a more rigorous and informative way.\nThe authors have responded to all the comments raised in the first round of review in a careful, systematic, and thorough manner. The revised manuscript shows significant improvements in terms of technical rigor, clarity of presentation, and overall completeness. The overall quality of the work now meets the standards for publication, and I recommend that the manuscript be accepted. Nevertheless, a few minor issues still require further clarification and refinement. I suggest that the authors provide more detailed explanations in the final version to further enhance the clarity and rigor of the manuscript.\n1. In the response to the reviewers’ comments, the authors provide a detailed timing breakdown for the DAC, ADC, and FPGA operations. Could the authors clarify how these timing statistics were obtained? Specifically, were they measured directly from hardware, or extracted from code level timing outputs or simulations?\nAuthors’ response: The timing statistics were measured using the FPGA Integrated Logic Analyzer (ILA). We probed relevant control and handshake signals and extracted cycle counts from captured waveforms. The reported values therefore correspond to cycle-accurate, on-chip measurements rather than simulation or code-level estimates. We have clarified this in the Method section.\n2. Regarding the description of emerging spintronic technologies for probabilistic computing, it would be beneficial to include references on related works, such as: Adv. Funct. Mater. 35, no. 43 (2025): 2501992. (https://doi.org/10.1002/adfm.202501992)\nAuthors’ response: We have added the suggested work.\nResponse to the fourth reviewer\nReviewer #4 (Remarks to the Author):\nThe authors have addressed all of the comments I had. The responses and the subsequent analysis are very comprehensive. I recommend acceptance of the manuscript. I believe the D-Wave implementation details are now sufficient to explain and motivate the quantum annealing sampling results. And the results that are described are consistent with the literature that addresses sampling using current D-Wave hardware.\nMy only small remaining comment is that the authors ran a NIST randomness test suite on the raw and post processed data. But that table and that data is not in the revised manuscript or the Supplementary information. I know this is not strictly the focus of the study, however, I do think it is interesting, and there could be folks interested in this technology for more security type applications, and so I would recommend putting this data somewhere (supplementary information, or even a data repository). I know that the techniques would need to be improved for truly secure TRNG applications, but I think this is a good experimental data to show somewhere. Or even just describing at least the reported p value results in written text would be sufficient. For example, quoting the min entropy in the SI document would be a very good reference point to state in my view; its a good information theoretic quantity that can be compared against other probabilistic generators.\nAuthors’ response: We thank the reviewer for this valuable suggestion. The NIST SP 800-22 test results, including detailed p-values and entropy analysis for the MTJ-based random number generator, were previously reported in Ref. [Ng et al., Phys. Rev. Appl. 19, 034077 (2023)]. As the randomness characterization has been comprehensively documented there, we cite that publication in the main text for readers interested in security-oriented performance metrics.\nThe authors have responded to all the comments raised in the first round of\nreview in a careful, systematic, and thorough manner. The revised manuscript shows\nsignificant improvements in terms of technical rigor, clarity of presentation, and\noverall completeness. The overall quality of the work now meets the standards for\npublication, and I recommend that the manuscript be accepted. Nevertheless, a few\nminor issues still require further clarification and refinement. I suggest that the\nauthors provide more detailed explanations in the final version to further enhance the\nclarity and rigor of the manuscript.\n1. In the response to the reviewers’ comments, the authors provide a detailed timing\nbreakdown for the DAC, ADC, and FPGA operations. Could the authors clarify how\nthese timing statistics were obtained? Specifically, were they measured directly from\nhardware, or extracted from code-level timing outputs or simulations?\n2. Regarding the description of emerging spintronic technologies for probabilistic computing,\nit would be beneficial to include references on related works, such as: Adv. Funct. Mater. 35,\nno. 43 (2025): 2501992. (https://doi.org/10.1002/adfm.202501992)"}, {"round": 2, "reviewer_comments": {"reviewer_2": "The author answered my question very well, and I recommend publishing the paper.", "reviewer_4": "The authors have addressed all of the comments I had. The responses and the subsequent analysis are very\ncomprehensive.\nI recommend acceptance of the manuscript.\nI believe the D-Wave implementation details are now sufficient to explain and motivate the quantum annealing sampling\nresults. And the results that are described are consistent with the literature that addresses sampling using current D-Wave\nhardware.\nMy only small remaining comment is that the authors ran a NIST randomness test suite on the raw and post processed data.\nBut that table and that data is not in the revised manuscript or the Supplementary information. I know this is not strictly the\nfocus of the study, however, I do think it is interesting, and there could be folks interested in this technology for more security\ntype applications, and so I would recommend putting this data somewhere (supplementary information, or even a data\nrepository). I know that the techniques would need to be improved for truly secure TRNG applications, but I think this is a\ngood experimental data to show somewhere. Or even just describing at least the reported p value results in written text\nwould be sufficient. For example, quoting the min entropy in the SI document would be a very good reference point to state\nin my view; its a good information theoretic quantity that can be compared against other probabilistic generators."}, "author_response": "The authors have responded to all the comments raised in the first round of\nreview in a careful, systematic, and thorough manner. The revised manuscript shows\nsignificant improvements in terms of technical rigor, clarity of presentation, and\noverall completeness. The overall quality of the work now meets the standards for\npublication, and I recommend that the manuscript be accepted. Nevertheless, a few\nminor issues still require further clarification and refinement. I suggest that the\nauthors provide more detailed explanations in the final version to further enhance the\nclarity and rigor of the manuscript.\n1. In the response to the reviewers’ comments, the authors provide a detailed timing\nbreakdown for the DAC, ADC, and FPGA operations. Could the authors clarify how\nthese timing statistics were obtained? Specifically, were they measured directly from\nhardware, or extracted from code-level timing outputs or simulations?\n2. Regarding the description of emerging spintronic technologies for probabilistic computing,\nit would be beneficial to include references on related works, such as: Adv. Funct. Mater. 35,\nno. 43 (2025): 2501992. (https://doi.org/10.1002/adfm.202501992)"}]

Decision letter