# Miniaturized time-correlated single-photon counting module for time-of-
flight non-line-of-sight imaging applications

Jie Wu  Hefei National Research Center for Physical Sciences at the Microscale
and School of Physical Sciences, University of Science and Technology of
China, Hefei 230026, China  CAS Center for Excellence in Quantum Information
and Quantum Physics, University of Science and Technology of China, Hefei
230026, China  Chao Yu<EMAIL_ADDRESS>Hefei National Research Center for
Physical Sciences at the Microscale and School of Physical Sciences,
University of Science and Technology of China, Hefei 230026, China  CAS Center
for Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei 230026, China  Jian-Wei Zeng  Hefei
National Research Center for Physical Sciences at the Microscale and School of
Physical Sciences, University of Science and Technology of China, Hefei
230026, China  CAS Center for Excellence in Quantum Information and Quantum
Physics, University of Science and Technology of China, Hefei 230026, China
Chen Dai  Hefei National Research Center for Physical Sciences at the
Microscale and School of Physical Sciences, University of Science and
Technology of China, Hefei 230026, China  CAS Center for Excellence in Quantum
Information and Quantum Physics, University of Science and Technology of
China, Hefei 230026, China  Feihu Xu  Hefei National Research Center for
Physical Sciences at the Microscale and School of Physical Sciences,
University of Science and Technology of China, Hefei 230026, China  CAS Center
for Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei 230026, China  Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China
Jun Zhang<EMAIL_ADDRESS>Hefei National Research Center for Physical
Sciences at the Microscale and School of Physical Sciences, University of
Science and Technology of China, Hefei 230026, China  CAS Center for
Excellence in Quantum Information and Quantum Physics, University of Science
and Technology of China, Hefei 230026, China  Hefei National Laboratory,
University of Science and Technology of China, Hefei 230088, China

###### Abstract

Single-photon time-of-flight (TOF) non-line-of-sight (NLOS) imaging enables
the high-resolution reconstruction of objects outside the field of view. The
compactness of TOF NLOS imaging systems, entailing the miniaturization of key
components within such systems is crucial for practical applications. Here, we
present a miniaturized four-channel time-correlated single-photon counting
module dedicated to TOF NLOS imaging applications. The module achieves
excellent performance with a 10 ps bin size and 27.4 ps minimum root-mean-
square time resolution. We present the results of TOF NLOS imaging experiment
using an InGaAs/InP single-photon detector and the time-correlated single-
photon counting module, and show that a 6.3 cm lateral resolution and 2.3 cm
depth resolution can be achieved under the conditions of 5 m imaging distance
and 1 ms pixel dwell time.

##  I Introduction

Non-line-of-sight (NLOS) imaging [1, 2, 3] is a novel technique for
reconstructing the details of obscured scenes from indirect light that has
been scattered multiple times. The ability to “see through walls” has
attracted extensive interest due to its potential applications in medical
detection, robotic surgery, earthquake relief, and autonomous driving. Over
the past decade, many studies have attempted to improve the performance of
NLOS imaging techniques in terms of resolution, distance, and imaging rate [4,
5, 6, 7, 8, 9]. Currently, the main implementation schemes for NLOS imaging
include time-of-flight (TOF) NLOS imaging [4], coherence-based NLOS imaging
[10], and intensity-based NLOS imaging [11]. In practice, TOF NLOS imaging
systems are widely used because of their robustness to the environment, low
cost, and high resolution.

Generally, the echo signal of a TOF NLOS imaging system is extremely weak
because of the multiple scattering. Since single-photon detector (SPD)
provides the ultimate sensitivity for photon detection, using SPD instead of
conventional photon detector can greatly improve the imaging sensitivity and
imaging distance of TOF NLOS imaging systems. Recently, TOF NLOS imaging
experiments using SPDs have been conducted for realistic applications [7, 8].
In practical scenarios, the single-photon TOF NLOS imaging system will be
integrated as a vehicle-mounted or handheld device, so the use of miniaturized
components is crucial. Among various SPD candidates, semiconductor SPDs based
on single-photon avalanche diodes have the advantages of small size, low cost,
and ease-of-use [12]. Recent studies have reported compact semiconductor SPDs
in different bands, including ultraviolet [13, 14], visible [15, 16], and
near-infrared [17, 18, 19], which could be highly suitable for single-photon
TOF NLOS imaging. Therefore, apart from the SPD module, it is also necessary
to miniaturize the remaining components, such as the time measurement module,
in single-photon TOF NLOS imaging systems.

In single-photon TOF NLOS imaging applications, the arrival time of an echo
signal is usually measured by a time-correlated single-photon counting (TCSPC)
module. Developing a highly integrated TCSPC module is crucial for practical
single-photon TOF NLOS imaging systems. There are several approaches to
achieve high-precision TCSPC. Conventionally, using a time-to-amplitude
convertor provides an analogue method of realizing picosecond time resolution
[20]. However, this method has inherent disadvantages, such as environmental
sensitivity, a short measurement range, and a relatively low count rate. In
contrast, digital methods, i.e., a time-to-digital converter (TDC), based on
the delay line architecture in field-programmable gate arrays (FPGAs) are
widely used due to the advantages of high resolution and flexible design [21,
22, 23, 24, 25]. Digital methods can also be implemented in application-
specific integrated circuits (ASICs) [26, 27]. Currently, high-performance
universal TDC ASIC chips are commercially available, which greatly helps the
development of an integrated TCSPC system.

Figure 1: (a) The principle of TCSPC technique. (b) The schematic diagram of
the whole application system with TCSPC, including a SPD, a TCSPC module and a
personal computer. (c) The photograph of the TCSPC testing board. (d) The
software interface for the TCSPC module.

In this paper, we present a miniaturized four-channel TCSPC module based on a
commercial TDC ASIC chip, i.e., TDC-GPX2 (Sciosense) [28]. An FPGA is used to
control the chip and handle the time measurement processing. A high-speed
universal serial bus (USB) 3.0 interface is used to transmit the time tag data
to a personal computer. A user-friendly software application allows time tag
data to be converted into a histogram. After characterizing the primary
parameters, including the timing resolution and non-linearity, we perform an
NLOS imaging experiment using the TCSPC module. The experimental results show
that the TCSPC module is suitable for the implementation of a compact single-
photon TOF NLOS imaging system.

##  II TCSPC MODULE DESIGN AND CHARACTERIZATION

The principle of the TCSPC technique is shown in Fig. 1 (a). A pulsed laser
illuminates the target with high repetition frequency. The echo photons are
collected and detected by an SPD, and the TCSPC electronics measure the
precise time difference between the photon arrival time and the reference
laser pulse. The time measurement data are then aggregated to form a
histogram, showing the intensity distribution of the echo photons at a given
time.

The whole application scheme using the TCSPC technique is depicted in Fig. 1
(b). The system consists of an SPD, TCSPC module and personal computer. In the
experiment, the SPD is implemented using a high-performance InGaAs/InP single-
photon avalanche diode [29]. To minimize the effect of stray photons, the SPD
operates in gated mode. The input gate clocks are synchronized to the laser
pulses, and the gate repetition frequency and gate width are set to 10 MHz and
10 ns, respectively. Following the standard calibration procedure [12], the
SPD exhibits excellent performance of 54% photon detection efficiency at 1550
nm, 2000 cps dark count rate, and 110 ps timing jitter. Moreover, the delay
time between the gates and the laser pulses can be tuned for each pixel during
the scanning process in the TOF NLOS imaging experiment.

The TCSPC module includes signal preprocessing, time measurement, logic
control, and data transmission functions, as shown in Fig. 1 (b). The module
supports input signals across four channels. In the TOF NLOS imaging
experiment, the four input signals are the start signal synchronized to the
laser pulses, SPD output signal, and two-dimensional pixel scan signals (pixel
$x$ and pixel $y$). The input signals are first discriminated, and then
converted into standard low voltage differential signaling (LVDS) with a fixed
pulse width of 10 ns using shapers. The LVDS signals are fed into TDC-GPX2 for
time measurement. An FPGA is used to initialize the TDC chip and read out the
time tag data, either in single-data-rate mode or in double-data-rate mode.
The FPGA also writes instructions to the chip registers via a serial
peripheral interface (SPI), allowing the operation modes of the TCSPC module
to be configured. The TCSPC module is implemented on an 8 cm $\times$ 12 cm
six-layer printed circuit board for testing purposes, in which the size of
time measurement component is $\sim$1 cm $\times$ 1 cm, as shown in Fig. 1
(c).

When the module is running, the arrival time of the input signal in each
channel is tagged based on the reference clocks generated by an external high-
accuracy crystal and a phase-locked loop (PLL) inside the FPGA. The bin size
is tuned by dividing the reference clocks with a configurable factor. The time
tag data consist of a coarse time package and a fine time package. The coarse
time is measured by counting the reference clock periods, while the fine time
is measured by counting the divided clock periods. During the experiment, the
frequency of the reference clock is set to 10 MHz. Considering the total
timing jitter of the NLOS imaging system is approximately 100 ps, the division
factor is initialized to 10000, corresponding to a 10 ps bin size. The length
of the fine time package is set to 14 bits, while the length of the coarse
time package is customized according to the repetition frequency of the NLOS
imaging system.

On the arrival of the input signals, the created time tag data are immediately
sent to the FPGA. To avoid data stacking, a memory architecture of first in
first out (FIFO) is designed inside the FPGA for data caching. In the TCSPC
module, input signals from channel 1 are defined as start signals, and input
signals from other channels are defined as stop signals. The time intervals
between the start and stop signals are calculated in real-time. The data are
then encoded with channel information, and transmitted to the personal
computer via the USB3.0 interface for data postprocessing. On the personal
computer, a user-friendly software is developed using PyCharm, as shown in
Fig. 1 (d). The software can flexibly configure various TCSPC parameters,
including the active channels, bin size, histogram depth, and measurement
time, and also periodically converts the time tag data to a histogram.

We then characterize the time resolution and non-linearity of the TCSPC
module. The time resolution is evaluated using the cable delay test technique
[30]. Pair signals under test are generated from the same clock signal with
different cable delays, and connected to channel 1 and channel 2 of the TCSPC
module. Figure. 2 (a) plots a typical time-interval histogram between paired
signals, which follows a Gaussian distribution. Given that the total number of
counts in the measurement is of order 1 million, the root-mean-square (RMS)
time resolution is 27.4 ps, as calculated by

$t_{RMS}=\sqrt{\frac{\sum_{i=1}^{n}{(t_{i}-t_{mean})}^{2}}{n-1}}.$ (1)

The measurements are repeated for various cable delays of up to 100 ns. As
shown in Fig. 2 (b), the RMS time resolution fluctuates from $\sim$27 ps to
$\sim$50 ps. In the NLOS imaging applications, the full width at half maximum
(FWHM) time resolution is usually used for calculating the spatial resolution.
By fitting the Gaussian distribution, the minimum FWHM time resolution of the
TCSPC module is calculated as 64.5 ps.

Figure 2: (a) Typical measured histogram of the time interval in the cable
delay test. The bin size is configured as 10 ps. (b) RMS time resolution
measurements by varying the cable delays.

The differential non-linearity (DNL) and integral non-linearity (INL) are
significant parameters for the TCSPC module. The DNL is the deviation of the
actual bin size from the ideal bin size, i.e., least significant bit (LSB).
The INL is calculated as the summation of the DNL from the first bin to the
current bin. We follow the method [31] to characterize the DNL and INL
parameters. The results are presented in Fig. 3 (a) and Fig. 3 (b),
respectively. The maximum values of the DNL and INL are 0.17 LSB and 8.49 LSB,
corresponding to 2 ps and 85 ps, respectively.

Figure 3: Characterization results of DNL (a) and INL (b). The bin size (LSB)
of the TCSPC module is configured as 10 ps.

##  III NLOS IMAGING EXPERIMENT

We then perform a single-photon TOF NLOS imaging experiment using the TCSPC
module. The experimental setup is shown in Fig. 4. A pulsed laser at 1550 nm
with 8 ps pulse width and 10 MHz repetition frequency is used to illuminate a
small patch of a visible diffuse wall, which is in the direct line of sight.
The light scatters off this patch to the hidden object, and then the object
reflects parts of the light back to the diffuse wall. In the confocal
receiving optical path, the reflected photons from the illuminated patch are
collected in a multimode fiber via a collimator and detected by the InGaAs/InP
SPD. Two two-axis galvanometers are used to scan an area of 0.8 m $\times$ 0.8
m on the diffuse wall with 64 $\times$ 64 pixels. The signal synchronized to
the laser pulses is connected to channel 1 of the TCSPC module as start
signal. The SPD output signal is connected to channel 2. The signals generated
by the galvanometer controller, which move the directions of pixel $x$ and
pixel $y$, are connected to channels 3 and 4, respectively. For each pixel,
the TCSPC module creates a TOF histogram. Using all scanning measurement data,
a three-dimensional image of the hidden target object is reconstructed
following the algorithm in Ref. [32].

Figure 4: Experimental setup for single-photon TOF NLOS imaging. Figure 5:
Reconstruction results for the four letter models of “USTC”

In the experiment, four white plastic letters “U”, “S”, “T”, and “C” with a
size of $\sim$ 0.5 m are used as hidden objects. The hidden object is $\sim$ 1
m away from visible diffuse wall. The distance between the visible diffuse
wall and the SPD is $\sim$ 5 m. As shown in Fig. 5, the four letters are
clearly reconstructed by the proposed scheme. The depth resolution
$\delta_{z}$ and lateral resolution $\delta_{x}$ of the imaging system are
estimated as 2.3 cm and 6.3 cm, respectively, as calculated by [7]

$\delta_{z}=c\sqrt{2Ln2(\frac{\sigma
cos(\phi)}{2c})^{2}+\frac{T_{SPD}^{2}+T_{laser}^{2}+T_{TCSPC}^{2}}{4}},$ (2)
$\delta_{x}=\delta_{z}\frac{\sqrt{w^{2}+z^{2}}}{w}.$ (3)

In Eqs. 2 and 3, $c$ represents the speed of light. $\sigma$ represents the
laser diameter, which is 2.67 cm. $\phi$ represents the angle between the
hidden object and the visible wall, which is 30∘. $T_{SPD}$, $T_{laser}$ and
$T_{TCSPC}$ represent the FWHM timing jitters of the SPD, laser, and the TCSPC
module, which are estimated as 110 ps, 8 ps and 64.5 ps, respectively. $z$ and
$w$ represent the distance from the hidden object to the visible wall and half
of the scanning width, which are 1 m and 0.4 m, respectively.

From Eqs. 2 and 3, the spatial resolution of the single-photon TOF NLOS
imaging system is primarily contributed by the optical parameters and the SPD
timing jitter. Even replacing the TCSPC module with a high-performance
commercial TCSPC instrument, e.g., the Time Tagger Ultra with 19 ps FWHM
timing jitter (Swabian Instrument) [33], the depth resolution and lateral
resolution can only be slightly improved to 2.2 cm and 5.8 cm, respectively,
while the size reaches 18 cm $\times$ 14 cm $\times$ 6 cm. From the
perspectives of integration and cost-effectiveness for applications, our TCSPC
module is highly suitable for practical single-photon TOF NLOS imaging
systems. In addition, the compactness and relatively simple hardware
implementation mean that the TCSPC module could be further integrated into the
InGaAs/InP SPD. Since the most functions of the TCSPC testing board, such as
power, signal preprocessing, logic control, and data transmission, are already
exist in the InGaAs/InP SPD, the additional size required for the integration
is estimated to be $\sim$4 cm2.

##  IV CONCLUSION

In summary, we have implemented a miniaturized TCSPC module dedicated to TOF
NLOS imaging applications with a commercially available TDC chip. The TCSPC
module exhibits excellent performance with a 10 ps bin size and 27.4 ps
minimum RMS time resolution. We have also demonstrated a single-photon TOF
NLOS imaging experiment using the TCSPC module, in which a 6.3 cm lateral
resolution and 2.3 cm depth resolution are achieved under the conditions of 5
m imaging distance and 1 ms pixel dwell time. Our work provides a practical
and cost-effective solution for the development of compact NLOS imaging
systems.

## acknowledgments

This work is supported by the Innovation Program for Quantum Science and
Technology under Grant Nos. 2021ZD0300804 and 2021ZD0300300, and the National
Natural Science Foundation of China (62175227).

## DATA AVAILABILITY

Data underlying the results presented in this paper are not publicly available
at this time but may be obtained from the authors upon reasonable request.

## REFERENCES

## References

  * [1] A. Kirmani, T. Hutchison, J. Davis, and R. Raskar, IEEE 12th International Conference on Computer Vision, 159-166 (2009).
  * [2] Y. Altmann, S. McLaughlin, M. J. Padgett, V. K. Goyal, A. O. Hero, and D. Faccio, Science 361, eaat2298 (2018).
  * [3] D. Faccio, A. Velten, and G. Wetzstein, Nat. Rev. Phys. 2, 318–327 (2020).
  * [4] A. Velten, T. Willwacher, O. Gupta, A. Veeraraghavan, M. G. Bawendi, and R. Raskar, Nat. Commun. 3, 745 (2012).
  * [5] A. Kirmani, D. Venkatraman, D. Shin, A. Colaço, F. N. C. Wong, J. H. Shapiro, and V. K. Goyal, Science 343, 58–61 (2014).
  * [6] M. Buttafava, J. Zeman, A. Tosi, K. Eliceiri, and A. Velten, Opt. Express 23, 1-7 (2015).
  * [7] C. Wu, J. Liu, X. Huang, Z.-P. Li, C. Yu, J.-T. Ye, J. Zhang, Q. Zhang, X. Dou, V. K. Goyal, F. Xu, and J.-W. Pan, Proceeding of the National Academy of Science 118, 1-7 (2021).
  * [8] B. Wang, M.-Y. Zheng, J.-J. Han, X. Huang, X.-P. Xie, F. Xu, Q. Zhang, and J.-W. Pan, Phys. Rev. Lett. 127, 053602 (2021).
  * [9] J. H. Nam, E. Brandt, S. Bauer, X. Liu, M. Renna, A. Tosi, E. Sifakis, and A. Velten. Nat. Commun. 12, 6526 (2021).
  * [10] C. A. Metzler, F. Heide, P. Rangarajan, M. M. Balaji, A. Viswanath, A. Veeraraghavan, and R. G. Baraniuk, Optica 7, 63-71 (2020).
  * [11] J. Klein, C. Peters, J. Martín, M. Laurenzis, and M. B. Hullin, Sci. Rep. 6, 32491 (2016).
  * [12] J. Zhang, M. A. Itzler, H. Zbinden, and J.-W. Pan, Light: Sci. Appl. 4, e286 (2015).
  * [13] C. Yu, T. Li, X.-S. Zhao, H. Lu, R. Zhang, F. Xu, J. Zhang, and J.-W. Pan, Rev. Sci. Instrum. 94, 033101 (2023).
  * [14] J.-T. Ye, C. Yu, W. Li, Z.-P. Li, H. Lu, R. Zhang, J. Zhang, F. Xu, and J.-W. Pan, Appl. Phys. Lett. 123, 024005 (2023).
  * [15] Y.-Q. Fang, K. Luo, X.-G. Gao, G.-Q. Huo, A. Zhong, P.-F. Liao, P. Pu, X.-H. Bao, Y.-A Chen, J. Zhang, and J.-W. Pan, Rev. Sci. Instrum. 91, 123106 (2020).
  * [16] F. Ceccarelli, G. Acconcia, A. Gulinatti, M. Ghioni, and I. Rech, IEEE Photonic. Tech. L. 31, 1 (2019).
  * [17] C. Yu, M. Shangguan, H. Xia, J. Zhang, X. Dou, and J.-W. Pan, Opt. Express 25, 14611–14620 (2017).
  * [18] W.-H. Jiang , X.-J. Gao, Y.-Q. Fang, J.-H. Liu, Y. Zhou, L.-Q. Jiang, W. Chen, G. Jin, J. Zhang, and J.-W. Pan, Rev. Sci. Instrum. 89, 123104 (2018).
  * [19] Q. Xu, C. Yu, W. Chen, J. Zhao, D. Cui, J. Zhang, and J.-W. Pan, arXiv:2310.17113 (2023).
  * [20] W. Becker, A. Bergmann, H. Wabnitz, D. Grosenick, and A. Liebert, European Conference on Biomedical Optics (2001).
  * [21] J.-Y. Wu, and Z.-H. Shi, IEEE Nuclear Science Symposium Conference Record, 3440-3446 (2008).
  * [22] Y. Wang, X. Zhou, Z. Song, J. Kuang, and Q. Cao, IEEE T. Nucl. Sci. 66, 2275-2281 (2019).
  * [23] D. Portaluppi, K. Pasquinelli, I. Cusini, and F. Zappa, IEEE T. Instrum. Meas. 71, 1-9 (2022).
  * [24] G. Cao, H. Xia, and N. Dong, Rev. Sci. Instrum. 89, 054707 (2018).
  * [25] X. Yu, H. Xia, W. Li, J. Zhang, and S. Chang, Rev. Sci. Instrum. 93, 085001 (2022).
  * [26] P. Dudek, S. Szczepanski, and J. V. Hatfield, IEEE J. Solid-St. Circ. 35, 240-247 (2000).
  * [27] C.-C. Chen, S.-H. Lin and C.-S. Hwang, IEEE T. Circuits-Ii. 61, 163-167 (2014).
  * [28] https://www.sciosense.com/products/time-to-digital-converters/tdc-gpx2-time-to-digital-converter/
  * [29] Y.-Q. Fang, W. Chen, T.-H. Ao, C. Liu, L. Wang, X.-J Gao, J. Zhang, and J.-W. Pan, Rev. Sci. Instrum. 91, 083102 (2020)
  * [30] J. Song, Q. An, and S. Liu, IEEE T. Nucl. Sci. 53, 236-241 (2006).
  * [31] J. Doernberg, H.-S. Lee, and D. A. Hodges, IEEE J. Solid-St. Circ. 19, 820 (1984).
  * [32] X. Liu, I. Guillen, M. L. Manna, J. H. Nam, S. A. Reza, T. H. Le, A. Jarabo, D. Gutierrez, and A. Velten, Nature 572, 620623 (2019).
  * [33] https://www.swabianinstruments.com/time-tagger/


# Digital Quantum Simulation of Reaction-Diffusion Systems on Lattice

Louie Hong Yao<EMAIL_ADDRESS>
###### Abstract

The quantum computer offers significant advantages in simulating physical
systems, particularly those with exponentially large state spaces, such as
quantum systems. Stochastic reaction-diffusion systems, characterized by their
stochastic nature, also exhibit exponential growth in the dimension of the
state space, posing challenges for simulation at a probability distribution
level. We explore the quantum simulation of stochastic reaction-diffusion
systems on a digital quantum computer, directly simulating the system at the
master equation level. Leveraging a spin representation of the system, we
employ Trotterization and probabilistic imaginary time evolution (PITE) to
simulate the probability distribution directly. We illustrate this approach
through four diverse examples, ranging from simple single-lattice site
generation-annihilation processes to a system featuring active-absorbing phase
transition.

##  I Introduction

Reaction-diffusion systems depict the dynamics of particles as they diffuse
through space and engage in local contact reactions. Despite their apparent
simplicity, these systems exhibit a wide array of fascinating phenomena,
including dynamical phase transitions[1, 2, 3, 4, 5, 6, 7, 8, 9], pattern
formation [10, 11, 12, 13, 14, 15], and self-organization [16, 17, 18].
Moreover, they serve as an effective framework for describing numerous
phenomena across diverse fields such as biology, ecology, chemistry,
economics, epidemiology, and game theory [19, 20, 21, 22, 23, 24, 25, 26, 27].

Traditionally, reaction-diffusion dynamics have been described using
interconnected nonlinear ordinary differential equations. To account for
spatially extended systems, this approach is extended to include partial
differential equations, incorporating species diffusion through simple
diffusion mechanisms. However, this traditional mean-field or mass-action
approach overlooks the inherent randomness and stochasticity in these
processes. These fluctuations, stemming from discrete individual numbers, and
spatio-temporal correlations are disregarded. Yet, these fluctuations and
correlations can lead to behaviors significantly different from those
predicted by mean-field theory [22].

The nonlinearities inherent in mean-field theory already pose challenges for
solving these systems, and the introduction of stochasticity further
complicates matters. Despite numerous attempts [28, 29, 30, 31], exact
solutions are rarely attainable due to the systems’ inherent stochasticity and
nonlinearity. Consequently, various analytical and numerical methods have been
employed to study them, including field theoretical approaches (see, for
example, [25, 11, 7, 3, 24]), agent-based Monte Carlo simulations (see, for
example, [26, 15, 32, 5]), and tensor network-related algorithms [33, 34, 35,
36, 37, 38]. However, most of these methods focus on understanding the system
at an averaged level by calculating the expectation values of observables of
interest. Although the master equation of the system is not difficult to write
down, the stochastic nature of the system leads to an exponentially growing
state space dimension. Consequently, conducting a probability distribution
level study within the framework of the master equation is challenging due to
constraints on computational power. Even for the simplest case with only a
single species of particles on lattice sites, the dimension of the state space
scales as $2^{N}$ with $N$ as the number of lattice sites. For $N>100$,
simulations become infeasible at all.

Quantum computers, leveraging qubits instead of classical bits, benefit from
the exponentially expanding Hilbert space as the number of qubits increases.
This unique feature has enabled them to demonstrate advantages in simulating
physical systems, especially those with exponentially large numbers of degrees
of freedom, such as quantum and stochastic systems. Various algorithms have
been devised to tackle these tasks. For instance, variational quantum
algorithms such as the variational quantum eigensolver [39, 40, 41] and
quantum approximate optimization algorithms [42, 43, 44] are designed to
compute the eigenspectrum of quantum systems through variational techniques.
Methods are also proposed to simulate quantum dynamics, in closed system [45,
46, 47, 48, 49, 50, 51, 52], open systems [53, 54, 55] and in imaginary time
for systems in thermal equilibrium [56, 57, 58, 59, 60, 61]. Algorithms have
also been developed to tackle partial differential equations [62, 63, 64, 65],
simulate stochastic systems [66, 67, 68, 69, 70, 71, 72], and explore fluid
dynamics [73, 74, 75]. The successful implementation of these algorithms on
real quantum devices (see, for example, [76, 77, 78, 79, 80]) addresses the
efficacy of quantum computing in simulating large-scale systems.

In this paper, we leverage the power of quantum computers to simulate
stochastic reaction-diffusion systems on lattices by directly simulating the
master equation governing the system. Our approach begins by mapping the
reaction-diffusion systems to spin models [29, 28], laying the groundwork for
direct simulation using qubits. The dynamics of the system are represented by
a spin ‘pseudo’-Hamiltonian or spin Liouvillian $H$, with the time evolution
governed by the operator $e^{-Ht}$. Since this ‘pseudo’-Hamiltonian is
generally non-anti-Hermitian, the time evolution is non-unitary. To simulate
this operator, we employ Trotterization methods. Additionally, we utilize the
probabilistic imaginary time evolution method proposed in [61] to implement
the non-unitary components within the Trotterization framework. To bridge the
gap in normalizations between quantum and stochastic states, we introduce
specialized pre- and post-processing techniques, eliminating the requirement
for exponential classical post-processing. We then put our methods to the test
across four diverse examples, spanning from simple single-site generation and
annihilation processes to a system demonstrating active-absorbing phase
transition. Through these examples, we showcase the effectiveness of our
simulation techniques in capturing various dynamic behaviors. While our
simulation results demonstrate promise in phases characterized by rapid system
relaxation, challenges emerge near critical points where slow relaxation
contributes to error accumulation, ultimately impacting result accuracy. In
our outlook, we propose potential solutions to address these limitations and
further enhance the robustness of the simulation results.

The article is structured as follows: In the subsequent section, we provide an
overview of reaction-diffusion dynamics and introduce spin representations of
these systems. Section III delves into our simulation methodologies, covering
Trotterization, Pauli gadget, probabilistic imaginary time evolution, and our
streamlined post-processing approach. We then proceed to apply these
simulation techniques to four distinct models in Section IV, where we compare
quantum simulation results with exact outcomes. Finally, we wrap up with a
concise summary and offer insights into future directions.

##  II A Brief Review of Reaction-Diffusion Dynamics and Spin Representations

Before we discuss the reaction-diffusion dynamics, we first review the master
equation description of the stochastic systems and present the pseudo-
Hamiltonian representation of the master equations. Given the probability
distribution $P$ of the system and the state space $\mathcal{H}$, most if not
all markovian stochastic dynamics can be written as

$\frac{\mathrm{d}}{\mathrm{d}t}P_{i}(t)=\sum_{j=1}^{D}W_{ij}P_{j}(t)-\sum_{j=1}^{D}W_{ji}P_{i}(t)$
(1)

where $D$ is the dimension of the state space and $i=1,2,3,...,D$ represent
the states in the state space $\mathcal{H}$. $W_{ij}$ are the transition rates
from state $j$ to state $i$.

Due to the linearity of the master equation, we can define
$|\Psi(t)\rangle=\sum_{i}P(i,t)|i\rangle$, where $\\{|i\rangle\\}$ are
orthogonal to each other. The master equation can then be written as a
‘pseudo’-Schrodinger equation:

$\frac{\mathrm{d}}{\mathrm{d}t}|\Psi(t)\rangle=-H|\Psi(t)\rangle,$ (2)

where the state $|\Psi(t)\rangle=\sum_{i}P_{i}(t)|i\rangle$ represents the
probability distribution of the system and $H$ is the pseudo-Hamiltonian (or
Liouvillian) with

$H|j\rangle=-\sum_{i}W_{ij}|i\rangle+\sum_{i}W_{ij}|j\rangle.$ (3)

The ‘pseudo’-Schrodinger equation Eq.(2) can be formally solved as

$|\Psi(t)\rangle=U(t)|\Psi(0)\rangle.$ (4)

with the time-evolution operator represented as $U(t)=e^{-Ht}$. Considering
observables represented by the operator $O$, the expectation values and two-
point correlation functions can be expressed as:

$\displaystyle\langle O(t)\rangle$ $\displaystyle=\langle
P|\hat{O}e^{-Ht}|\Psi(0)\rangle,$ (5) $\displaystyle\langle
O(t_{2})O(t_{1})\rangle$ $\displaystyle=\langle
P|\hat{O}e^{-H(t_{2}-t_{1})}\hat{O}e^{-Ht_{1}}|\Psi(0)\rangle$

where $\langle P|=(1,\ 1,\ ...,1)$ is a projection state. The normalization
condition of the probability is then given by $\langle P|\Psi\rangle=1$.

The reaction-diffusion dynamics involve the diffusion of particles and their
reactions when they come into contact. A simple example is pair annihilation
dynamics, where particles of the same type move around and may vanish when
they collide. In this paper, we’re focusing on systems arranged on lattices
with just one species of particles allowed. Our simulation methods can be
extended to systems with multiple species, either by using qudits or by using
multiple qubits to represent one lattice site. We also enforce a rule that
each lattice site can only hold one particle, which is typical in Monte Carlo
simulations for these systems. Each lattice site then has two possible states:
occupied (denoted as $|\uparrow\rangle$) or empty (denoted as
$|\downarrow\rangle$), making the state space for a single site
$\mathcal{H}_{2}=\mathbb{R}^{2}$. Extending this concept to a lattice
comprising $N$ sites yields the comprehensive state space
$\mathcal{H}=\mathcal{H}_{2}^{\otimes N}$. Particle generation and
annihilation are succinctly captured by the operators
$\sigma^{\pm}=(\sigma^{x}\pm i\sigma^{y})/2$, following established
conventions of quantum mechanics. Hence, particle interactions can be depicted
by utilizing these operators. A fundamental instance arises in the
annihilation of a single particle with a rate $\lambda$. The Hamiltonian
governing this single-particle reaction is depicted as
$\lambda\sum_{i}(n_{i}-\sigma_{i}^{-})$, where
$n_{i}=\sigma_{i}^{+}\sigma_{i}^{-}=(\sigma_{i}^{z}+1)/2$ represents the
particle number operator, and the index $i$ denotes the lattice sites. Other
reactions can also be expressed in this way. We outline the Hamiltonians
corresponding to common reactions in Table 1. In cases involving multiple
reactions, the Hamiltonian for the entire system is derived by summing the
Hamiltonians of each reaction. For example, the system with particle diffusion
and pair annihilation can be represented by

$\displaystyle H=$ $\displaystyle D\sum_{\langle
i,j\rangle}\Bigg{[}(\sigma_{i}^{+}-\sigma_{j}^{+})(\sigma_{i}^{-}-\sigma_{j}^{-})-2n_{i}n_{j}\Bigg{]}$
(6) $\displaystyle+\lambda\sum_{\langle
i,j\rangle}(n_{i}n_{j}-\sigma_{i}^{-}\sigma_{j}^{-}),$

where $\langle i,j\rangle$ denotes nearest neighbors. Given that the three
Pauli matrices $\sigma^{x},\ \sigma^{y},\ \sigma^{z}$ and the 2-dimensional
identity matrix constitute a complete basis for $SL(2,\mathbb{C})$, the
Hamiltonian can be represented as linear combinations of tensor products of
the Pauli basis $\\{1,\sigma^{x},\sigma^{y},\sigma^{z}\\}$. The expression of
the Hamiltonian in the Pauli basis is also included, as exponentials of Pauli
matrices are simpler to implement on a digital quantum computer.

| Pseudo-Hamiltonians
---|---
Reactions | Generation-Annihilation basis | Pauli basis
Free particles Hopping | $\sum_{\langle i,j\rangle}(\sigma_{i}^{+}-\sigma_{j}^{+})(\sigma_{i}^{-}-\sigma_{j}^{-})-2n_{i}n_{j}$ | $-\frac{1}{2}\sum_{\langle i,j\rangle}\Big{(}\sigma_{i}^{x}\sigma_{j}^{x}+\sigma_{i}^{y}\sigma_{j}^{y}+\sigma_{i}^{z}\sigma_{j}^{z}-1\Big{)}$
Pair Annihilation ($2A\rightarrow 0$) | $\sum_{\langle i,j\rangle}(n_{i}n_{j}-\sigma_{i}^{-}\sigma_{j}^{-})$ | $\begin{aligned} \frac{1}{4}\sum_{\langle i,j\rangle}\Big{(}-\sigma_{i}^{x}\sigma_{j}^{x}+\sigma_{i}^{y}\sigma_{j}^{y}+\sigma_{i}^{z}\sigma_{j}^{z}+i\sigma_{i}^{x}\sigma_{j}^{y}+i\sigma_{i}^{y}\sigma_{j}^{x}\Big{)}\\\ +\frac{1}{2}\sum_{i}\left(\sigma_{i}^{z}+\frac{1}{2}\right)\end{aligned}$
Pair Coagulation ($2A\rightarrow A$) | $\sum_{\langle i,j\rangle}n_{i}n_{j}-\frac{1}{2}n_{i}\sigma_{j}^{-}-\frac{1}{2}n_{j}\sigma_{i}^{-}$ | $\begin{aligned} \frac{1}{8}\sum_{\langle i,j\rangle}\Big{(}2\sigma_{i}^{z}\sigma_{j}^{z}-\sigma_{i}^{z}\sigma_{j}^{x}-\sigma_{i}^{x}\sigma_{j}^{z}+i\sigma_{i}^{z}\sigma_{j}^{y}+i\sigma_{i}^{y}\sigma_{j}^{z}\Big{)}\\\ +\frac{1}{4}\sum_{i}\left(2\sigma_{i}^{z}-\sigma_{i}^{x}+i\sigma_{i}^{y}+\frac{1}{4}\right)\end{aligned}$
Spontaneous decay ($A\rightarrow 0$) | $\sum_{i}n_{i}-\sigma_{i}^{-}$ | $\frac{1}{2}\sum_{i}\Big{(}-\sigma_{i}^{x}+i\sigma_{i}^{y}+\sigma_{i}^{z}+1\Big{)}$
Spontaneous generation ($0\rightarrow A$) | $\sum_{i}(1-n_{i})-\sigma_{i}^{+}$ | $\frac{1}{2}\sum_{i}\Big{(}1-\sigma_{i}^{z}-\sigma_{i}^{x}-i\sigma_{i}^{y}\Big{)}$
Branching ($A\rightarrow 2A$) | $\sum_{\langle i,j\rangle}n_{i}+n_{j}-2n_{i}n_{j}-n_{i}\sigma_{j}^{+}-n_{j}\sigma_{i}^{+}$ | $\begin{aligned} -\frac{1}{4}\sum_{\langle i,j\rangle}\Big{(}2\sigma_{i}^{z}\sigma_{j}^{z}+\sigma_{i}^{z}\sigma_{j}^{x}+\sigma_{i}^{x}\sigma_{j}^{z}+i\sigma_{i}^{z}\sigma_{j}^{y}+i\sigma_{i}^{y}\sigma_{j}^{z}\Big{)}\\\ -\frac{1}{2}\sum_{i}\left(\sigma_{i}^{x}+i\sigma_{i}^{y}-1\right)\end{aligned}$
Table 1: The Pseudo-Hamiltonians for common reactions involving particle
species $A$ are detailed above. Expressions are provided both in terms of
$\\{\sigma^{+},\sigma^{-},n\\}$ and
$\\{1,\sigma^{x},\sigma^{y},\sigma^{z}\\}$. As we only consider contact
reactions and nearest-neighbor hopping, $\langle i,j\rangle$ represents
nearest-neighbor pairs.

##  III Methods

In this section, we delve into the simulation methods employed to simulate the
reaction-diffusion system on a digital quantum computer. Additionally, we
elaborate on several implementation details pivotal to our simulation process.

###  III.1 Trotterization

Quantum simulation of dynamics, whether in real-time, imaginary time, or
involving general open quantum dynamics, requires implementing operators of
the form $e^{-Ht}$, where $H$ represents the system’s Hamiltonian,
Liouvillian, or Lindbladian. Even if $H$ consists solely of local terms, i.e.,
$H=\sum_{i}H_{i}$ where $H_{i}$s are local, the operator $e^{-Ht}$ can still
pose significant nonlocality challenges for implementation. Various simulation
methods have been proposed to tackle this challenge, including the Lie-
Trotter-Suzuki decomposition [48, 49], linear combination of unitaries (LCU)
[50, 52, 81], and quantum signal processing (QSP) [51]. In this paper, we
choose trotterization due to its intuitive nature, ease of implementation, and
demonstrated good empirical performance [82].

Given the challenge of handling the operator $e^{-Ht}$, our aim is to
decompose it into a product $\prod_{i}e^{-H_{i}t}$, where $e^{-H_{i}t}$ terms
are more manageable. In many physical systems, including the reaction-
diffusion systems under discussion, $H$ can naturally be expressed as
$H=\sum_{i}H_{i}$, where $H_{i}$s represent local terms and are easy to
implement.

However, due to the non-commutativity between $H_{i}$s, $e^{-\sum_{i}H_{i}t}$
isn’t strictly equal to $\prod_{i}e^{-H_{i}t}$. An effective solution is
offered by the Lie product formula, where given two arbitrary matrices $A$ and
$B$, we have:

$e^{A+B}=\lim_{n\rightarrow\infty}(e^{\frac{A}{n}}e^{\frac{B}{n}})^{n}.$ (7)

Thus, we can divide $e^{-Ht}$ into small trotter steps $\Delta t$ and express
it as:

$e^{-Ht}=\lim_{\Delta t\rightarrow 0}\left(\prod_{i}e^{-H_{i}\Delta
t}\right)^{\frac{t}{\Delta t}}.$ (8)

This concept of trotterization has been utilized by physicists for some time,
including in the derivation of coherent state path integrals [83, 25].
However, in the implementation on a quantum device, we cannot take the limit
$\Delta t\rightarrow\infty$. It has been shown that, by the
Baker–Campbell–Hausdorff formula, to first order [48]:

$e^{-\sum_{i}H_{i}t}=\left(\prod_{i}e^{-H_{i}\Delta t}\right)^{\frac{t}{\Delta
t}}+\mathcal{O}\left(\frac{t^{2}}{n}\right),$ (9)

where $n=t/\Delta t$. Thus, by choosing small enough $\Delta t$, we can
approximate the real dynamics $e^{-Ht}$ very effectively. A more detailed
error analysis of this method can be found in [84]. Throughout this paper, we
will only consider the first-order approximation.

Other methods exist for decomposing $e^{-Ht}$ exactly, such as the Cartan
decomposition [85]. However, these methods require extensive classical
preprocessing and we will not utilize them in this work.

###  III.2 Pauli Gadgets

The Hamiltonian $H$ of spin systems and reaction-diffusion systems with site-
restriction can be expressed as a linear combination of tensor products of
Pauli matrices. This implies that the $H_{i}$s are tensor products of Pauli
matrices. Implementing these operations on quantum computers necessitates
expressing these gates using universal gate sets, typically Clifford gates
with a single non-Clifford gate. One straightforward method is through the use
of Pauli gadgets, which initially transform the operation to the $z$-basis via
rotations in single-qubit Hilbert space. Then, information is propagated to a
single target qubit, where a special single-qubit rotation is applied. For
example, consider the operation
$e^{-i\sigma^{x}\otimes\sigma^{y}\otimes\sigma^{z}t}$. We can express it as:

$e^{-i\sigma^{x}\otimes\sigma^{y}t}=(H\otimes S)(H\otimes
I)e^{-i\sigma^{z}\otimes\sigma^{z}t}(H\otimes H)(S^{\dagger}\otimes I),$ (10)

where

$H=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\ 1&-1\end{pmatrix},\qquad
S=\begin{pmatrix}1&0\\\ 0&i\end{pmatrix}$ (11)

are the Hadamard gate and the $S$ gate respectively. With this transformation,
we then only need to implement
$e^{-i\sigma^{z}\otimes\sigma^{z}\otimes\sigma^{z}t}$.

Operations $T\left((\sigma^{z})^{\otimes n}\right)$, composed solely of tensor
products of Pauli Z matrices, can be expressed as
$\sum_{i}T(p_{i})|i\rangle\langle i|$, where $i$ represents binary strings.
Here, $p_{i}$ equals $1$ when the binary string $i$ contains an even number of
$1$s, and $p_{i}$ equals $-1$ when the binary string $i$ contains an odd
number of $1$s. We then utilize a CNOT gate to propagate the information to a
single qubit, implementing an $R_{z}$ gate based on the state (0 or 1) of that
qubit. For the example $e^{-i\sigma^{x}\otimes\sigma^{y}\otimes\sigma^{z}t}$,
we propagate to the third qubit. The circuit implementation is depicted in
Fig. 1. As plotted in the circuit, we initially employ the $H$ and $S$ gates
to transform to the $z$-basis. Subsequently, we utilize CNOT gates to transmit
information to the third qubit and apply an $R_{z}$ gate, i.e.

$R_{z}(\theta)=e^{-i\frac{\theta}{2}\sigma^{z}}=\begin{pmatrix}e^{-i\frac{\theta}{2}}&0\\\
0&e^{i\frac{\theta}{2}}\end{pmatrix}.$ (12)

Following these steps, we revert the system back via inverse transformations.
It’s essential to note that quantum circuits are applied from left to right,
contrary to the usual Dirac notation representation.

In this entire implementation, the Hadamard gate, the $S$ gate, and the CNOT
gate are all Clifford gates. The entire procedure only necessitates a single
non-Clifford gate, $R_{z}$. Since Clifford gates are simpler to implement on
most current quantum devices, the Pauli gadgets method significantly
simplifies the implementation process at the hardware level. Further
optimizations of the circuit implementations via Pauli gadgets has also been
done and can be found in, for example [86, 87].

Employing Pauli gadgets offers another advantage: the reduction in the number
of ancilla qubits required for non-unitary operations [61]. During the Pauli
gadgets process, information is propagated to a single qubit, necessitating
only a single-qubit non-unitary operation. It is well-known that we can purify
a system with Hilbert space $\mathcal{H}$ in a tensor product Hilbert space
$\mathcal{H}\otimes\mathcal{H}$. Consequently, only a single ancilla qubit is
needed to implement the non-unitary operation of the form $e^{-Ht}$ when $H$
is a tensor product of Pauli matrices.

[node distance=1.5pt] &H 1 1 H
S^† H1 1 H S
R_Z(2t)

Figure 1: Implementing the time evolution operator
$e^{-i\sigma^{x}\otimes\sigma^{y}\otimes\sigma^{z}t}$ through Pauli gadgets
involves a sequence of steps. Initially, we rotate the system to the Pauli $z$
basis. Subsequently, we propagate information to the third qubit and apply an
$R_{z}$ gate. Finally, we perform the inverse operation to return the system
to its original basis.

###  III.3 Probabilistic Imaginary Time Evolution

Since the pseudo-Hamiltonian $H$ is not generally anti-Hermitian, simulating
the time evolution $e^{-Ht}$ requires implementing non-unitary operations.
According to Trotterization to first order, $e^{-Ht}$ can be decomposed into
$(\prod_{i}e^{-H_{i}t/N})^{N}$, where $H_{i}$ can be either Hermitian or anti-
Hermitian. For anti-Hermitian $H_{i}$, $e^{-H_{i}t/N}$ is unitary and thus
straightforward to implement. Conversely, for Hermitian $H_{i}$, implementing
$e^{-H_{i}t/N}$ is referred to as imaginary time evolution, as it is
equivalent to $e^{-iH_{i}\tau}$ with $\tau=-it$. The time $t$ here can also be
understood as the effective temperature $\beta$ as the imaginary time
evolution corresponds to the equilibrium partition function of Gibbs
distribution.

Several methods have been proposed for simulating imaginary time evolution on
quantum computers. Variational approaches, such as variational imaginary time
evolution (VITE) [56, 57] and quantum imaginary time evolution (QITE) [58, 59,
60], offer a hybrid approach by optimizing a variational ansatz to approximate
the target state as closely as possible. Another approach, the probabilistic
approach [88, 89, 90, 91, 92, 93, 61], which appears more natural to
physicists, involves introducing ancilla qubits or environments to the system
and then applying post-selection while measuring out the ancilla qubits. In
this work, we adapt the methods proposed in [61], which leverage the power of
the Pauli gadgets introduced in Section III.2 to reduce the number of ancilla
qubits used in the system.

With the Pauli gadgets, the Trotterized real-time evolution operator can be
realized by single-qubit gates $e^{-i\alpha\sigma^{z}}$ with single-qubit
rotations and two-qubit control gates, as shown in Fig.1. Similarly, for the
imaginary time evolution operator, we can map it to a single-qubit non-unitary
operation $e^{-\alpha\sigma^{z}}$. Thus, we only need to implement the single-
qubit non-unitary operation of the form $e^{-\alpha\sigma^{z}}$. To do so, we
can introduce a single ancilla qubit and implement a controlled-RX gate
($CR_{x}$). The $R_{x}$ gate is given by

$R_{x}(\theta)=e^{-i\frac{\theta}{2}\sigma^{x}}=\begin{pmatrix}\cos\frac{\theta}{2}&-i\sin\frac{\theta}{2}\\\
-i\sin\frac{\theta}{2}&\cos\frac{\theta}{2}\end{pmatrix}.$ (13)

Given a state $a|0\rangle+b|1\rangle$ and an ancilla qubit at $|0\rangle$, the
$CR_{x}$ gate transforms the two-qubit state to

$CR_{x}\left(a|00\rangle+b|10\rangle\right)=a|00\rangle+b\cos\frac{\theta}{2}|10\rangle-
ib\sin\frac{\theta}{2}|11\rangle,$ (14)

where the second qubit is the ancilla qubit. By measuring the second qubit and
post-selecting it to be $0$, we obtain the final state
$a|0\rangle+b\cos\frac{\theta}{2}|1\rangle$ up to a normalization factor.
Since the non-unitary operation $e^{-\alpha\sigma^{z}}$ yields

$e^{-\alpha\sigma^{z}}(a|0\rangle+b|1\rangle)\propto
a|0\rangle+e^{2\alpha}b|1\rangle,$ (15)

when $\alpha<0$, we can implement the $CR_{x}$ gate with
$\theta=2\arccos\left(\exp(-2|\alpha|)\right)$ and post-select the ancilla
qubit to be $0$. This procedure is depicted by the circuit shown in Fig.2(a).
When $\alpha>0$, $e^{2\alpha}$ is larger than $1$ and cannot be mapped to a
cosine. Instead, we can use the Pauli $X$ gate to bring $\cos\frac{\theta}{2}$
to the state $|0\rangle$, i.e.,

$\displaystyle(\sigma^{x}\otimes 1)CR_{x}$ $\displaystyle(\sigma^{x}\otimes
1)\left(a|00\rangle+b|10\rangle\right)$ (16)
$\displaystyle=a\cos\frac{\theta}{2}|00\rangle+b|10\rangle-
ia\sin\frac{\theta}{2}|01\rangle.$

By selecting $\theta=2\arccos\left(\exp(-2|\alpha|)\right)$, we obtain the
correct outcome. The circuit for $\alpha>0$ is shown in Fig. 2(b). With
Trotterization and Pauli gadgets combined, we are now equipped to implement
any time evolutions with the ’Pseudo’-Hamiltonian listed in Table 1.

[node distance=1.5pt,wire types=q,n,column sep=10pt, row sep=20pt,between
origins] & 1
$\ket{0}$ q R_x(θ) 00n

(a)

[node distance=1.5pt,wire types=q,n,column sep=10pt, row sep=20pt,between
origins] & X1 X
$\ket{0}$ q R_x(θ) 00n

(b)

Figure 2: The quantum circuit implementation of non-unitary operation
$e^{-\alpha\sigma^{z}}$ for (a) $\alpha<0$ and (b) $\alpha>0$. The angle
$\theta$ is given by $\theta=2\arccos\left(\exp(-2|\alpha|)\right)$.

###  III.4 Encoding and Post-processing

While the pseudo-Schrödinger Eq.(2) can be represented by combined unitary and
non-unitary gates in a quantum circuit, there are still distinctions between
this formulation and the quantum-mechanical Schrödinger’s equation. For
clarity, we denote the quantum states in lowercase and the probability state
in uppercase.

One notable difference lies in the normalization of states. In quantum
mechanics, state vectors $|\psi\rangle$ are $L_{2}$ functions with the
$2$-norm normalized to $1$, meaning $\langle\psi|\psi\rangle=1$. However, in
the study of reaction-diffusion systems, states are normalized as $\langle
P|\Psi\rangle=1$, or $\sum_{i}P_{i}=1$, in accordance with probability
conservation. Consequently, necessary encoding and decoding operations must be
performed.

For encoding, we can simply rescale the initial state as follows:

$|\Psi\rangle=\sum_{i}P_{i}|i\rangle\quad\Longrightarrow\quad|\psi\rangle=\frac{1}{\sqrt{\sum_{i}P_{i}^{2}}}\sum_{i}P_{i}|i\rangle.$
(17)

After the time-evolution, the resulting state from the quantum circuit is
expressed as:

$|\psi_{t}\rangle=\frac{1}{\sqrt{\langle\psi|U^{\dagger}_{t}U_{t}|\psi\rangle}}U_{t}|\psi\rangle,$
(18)

where $U_{t}=e^{-Ht}$ represents the time evolution operator. To ensure
correct normalization, we need to compute $\langle P|\psi_{t}\rangle$. The
projection state $\langle P|$ can be represented as:

$\langle P|=\sqrt{2^{n}}\langle\mathbf{0}|H^{\otimes n},$ (19)

assuming there are $n$ qubits in the system. Here $\langle\mathbf{0}|$
signifies the all-zero state $\langle 0|^{\otimes n}$. As a result, the
normalization factor $\langle P|\psi_{t}\rangle$ is given by:

$\langle P|\psi_{t}\rangle=\sqrt{2^{n}}\langle\mathbf{0}|H^{\otimes
n}|\psi_{t}\rangle.$ (20)

By measuring $\langle\mathbf{0}|H^{\otimes n}|\psi_{t}\rangle$, we obtain the
probability amplitude of observing all qubits in the state $H^{\otimes
n}|\psi_{t}\rangle$ to be in the state $|0\rangle$. Essentially, this quantity
gives us the amplitude for observing all zeroes when measuring the state
$H^{\otimes n}|\psi_{t}\rangle$. Since this measurement outcome can be
directly obtained from the final measurement of the quantum circuit, it
eliminates the need for any normalization calculations on a classical
computer, which could otherwise be computationally expensive due to its
exponential nature.

Since both the normalization factor $\langle P|\psi_{t}\rangle$ and all
components of the state $|\psi_{t}\rangle$ are real-valued, they can be
obtained from measurements of probabilities alone. Consequently, there’s no
need for additional tomography or complex calculations; everything can be
directly inferred from the experimental outcomes.

##  IV Models and Results

In this section, we apply the methods introduced in Section III to four
distinct examples on a one-dimensional lattice. These examples span from a
straightforward scenario involving single particle generation and annihilation
to a model showcasing an active-absorbing phase transition. Given the
limitations of current quantum devices, we perform classical simulations of
the process and compare the outcomes with the exact results.

###  IV.1 Single Particle Generation and Annihilation

We begin with a simple example involving only a single lattice site. Here, we
consider single particle generation and annihilation reactions on the site
with rates $\nu$ and $\lambda$, respectively. The ‘pseudo’-Hamiltonian for
this system is given by

$H=\frac{1}{2}\left[-(\lambda+\nu)\sigma^{x}+(\lambda-\nu)\sigma^{z}+i(\lambda-\nu)\sigma^{y}\right],$
(21)

where we have omitted the lower indices for the lattice site. Additionally, we
have omitted the constant terms since they only contribute to rescaling in the
time-evolution. Neglecting these constants should not affect the results, as
we will need to rescale the outcome from the quantum simulation at the end of
the process. We will automatically disregard these constant terms in the
subsequent calculations unless explicitly stated otherwise. In this system,
particle generation competes with annihilation. In the limit $\nu=0$, the
stationary state will be an empty site. Conversely, in the limit where
$\lambda=0$, the particle density of the stationary state will be $1$. For
values between these two limits, we will observe a finite particle density
$0<\langle n\rangle<1$ in the stationary state.

We simulate the process and depict the time evolution of the particle density
$\langle n\rangle$ starting from an occupied state, as shown in Fig.(3). In
the simulation, the ’pseudo’-Hamiltonian is trotterized into terms
$e^{(\lambda+\nu)\sigma^{x}\Delta t}$, $e^{-(\lambda-\nu)\sigma^{z}\Delta t}$,
and $e^{-i(\lambda-\nu)\sigma^{y}\Delta t}$, where the first two terms are
non-unitary and the last term is unitary and can be implemented directly. For
simplicity, we choose $\lambda=1$ in the simulation and compare the results
for different values of $\nu$.

In Fig. 3(a), we plot the particle density $\langle n\rangle$ versus time $t$
for different values of $\nu$ with a trotter step size $\Delta t=1/20$. The
solid lines represent the exact results, while the scattered labels represent
the results from quantum simulations. We observe that for intermediate values
of $\nu$, the process is very well captured by the simulation. However, for
both large and small values of $\nu$, there are visible errors, which arise
from the fact that either $\lambda+\nu$ or $\lambda-\nu$ is large, leading to
increased trotter errors.

We further simulate the cases of large and small $\nu$ values with a smaller
trotter step size $\Delta t=1/50$ in Fig.3(b). As we decrease the trotter step
size, the errors decrease, and the results better agree with the exact
results.

(a)

(b)

Figure 3: The dynamics of particle density in the single-site generation and
annihilation system with fixed generation rate $\lambda=1$ and varying
annihilation rate $\nu$s: In (a), the trotter step size is $\Delta t=1/20$,
while in (b), it is $\Delta t=1/80$. The solid lines represent the exact
results, while the markers indicate the simulated results.

###  IV.2 Free Particles Hopping

The next example we consider involves a free particle hopping on a line. Given
the diffusion constant $D$, the ‘pseudo’-Hamiltonian is written as

$H=-\frac{1}{2}D\sum_{\langle
i,j\rangle}\left(\sigma_{i}^{x}\sigma_{j}^{x}+\sigma_{i}^{y}\sigma_{j}^{y}+\sigma_{i}^{z}\sigma_{j}^{z}\right).$
(22)

This ‘pseudo’-Hamiltonian coincides with the Hamiltonian of the Heisenberg
model; however, the time evolution is not unitary. We benchmark the simulation
on a system with $4$ lattice sites, and the results are plotted in Fig.4.
Since particle hopping conserves the total particle number in the system, we
instead plot the dynamics of the probability of the relevant states. The
trotter step size is chosen as $\Delta t=1/20$.

In Fig.4(a), we simulate the system starting with the state
$\bigstar\bigcirc\bigcirc\bigcirc$ with a diffusion constant $D=1$. Here,
$\bigstar$ indicates the site is occupied, and $\bigcirc$ indicates the site
is empty. Since there is only one single particle in the system, only four
states are relevant in the dynamics, as shown in the plot. As expected, we
observe that the probability of all four states converges to $1/4$.
Additionally, due to the periodic boundary condition we choose, the states
labeled with a triangle and a plus sign are symmetric to each other.
Consequently, their probabilities are equal, leading to an overlap in the
plot. In Fig.4(b), we present the results from a different initial state:
$2/3$ probability of staying in $\bigstar\bigcirc\bigcirc\bigcirc$ and $1/3$
probability of staying in $\bigcirc\bigstar\bigcirc\bigcirc$. Here, we employ
a different diffusion constant $D=0.6$ in the simulation. As anticipated, we
observe that the probability of all four relevant states converges to $1/4$,
resulting in a uniform distribution. This uniform distribution is in line with
the expected behavior of a system undergoing free particle hopping, where the
particle density tends to spread out evenly across the lattice over time. From
the plots, we can see that at a trotter step size $\Delta t=1/20$, the
simulation recovers the exact result very well at the probability density
level. This indicates that the trotterization method, combined with the Pauli
gadgets, effectively captures the non-unitary evolution described by the
’pseudo’-Hamiltonian.

(a)

(b)

Figure 4: The probabilities associated with different states in the pure
hopping dynamics on four lattice sites: In (a), the initial state is
$|\bigstar\bigcirc\bigcirc\bigcirc\rangle$, and the diffusion constant is
$D=1$. In (b), the initial state is
$2/3|\bigstar\bigcirc\bigcirc\bigcirc\rangle+1/3|\bigcirc\bigstar\bigcirc\bigcirc\rangle$,
and the diffusion constant is $D=0.6$. Here $\bigstar$ represents the occupied
site and $\bigcirc$ represents the empty site.

###  IV.3 Pair Annihilation

In the third example, we consider the pair annihilation reaction
$2A\rightarrow 0$ with a reaction rate $\nu$, alongside particles being able
to freely hop on the lattice with a diffusion constant $D$. The
‘pseudo’-Hamiltonian governing the dynamics is given by:

$\displaystyle H=-\frac{1}{4}\sum_{\langle
i,j\rangle}\Big{[}(2D+\nu)\sigma_{i}^{x}\sigma_{j}^{x}+(2D-\nu)\sigma_{i}^{y}\sigma_{j}^{y}$
(23)
$\displaystyle+(2D-\nu)\sigma_{i}^{z}\sigma_{j}^{z}-i\nu\sigma_{i}^{x}\sigma_{j}^{y}-i\nu\sigma_{i}^{y}\sigma_{j}^{x}\Big{]}$
$\displaystyle+\frac{\nu}{2}\sum_{i}\sigma_{i}^{z}.$

For simplicity, we set the diffusion constant $D=1$ in the simulation. Due to
the requirement of the reaction that two particles are needed to annihilate,
the total number of particles initially in the system affects the final state.
We study two different cases: a fully-occupied six-site lattice and a fully-
occupied seven-site lattice, as plotted in Fig.5.

In Fig.5(a), we display the total particle number $\langle n\rangle$ in the
six-lattice site with a fully occupied initial condition. For different values
of $\nu$, the system converges to an empty state at different rates. Fig.5(b)
illustrates the time evolution of the total particle number in a seven-lattice
site system with a fully occupied initial condition. Unlike the case with $6$
particles initially, the system converges to a final state with $1$ particle
remaining, as it cannot find a partner to annihilate itself. In the seven-site
system, the trotter step size is chosen as $\Delta t=1/50$, and in the six-
site system, the trotter step size is chosen as $\Delta t=1/200$.
Surprisingly, to achieve the same level of accuracy, we require a smaller
trotter step size in the six-site system compared to the seven-site system.
The challenge arises from the nature of the systems being simulated. In the
case of the six-site system, it converges to a completely empty state, whereas
the seven-site system converges to a superposition of different single-
particle states. In such scenarios, the stationary state, representing the
final outcome of the system, should correspond to one of the right eigenstates
of the ’pseudo’-Hamiltonian. However, during the Trotterization process, the
empty state might not be an eigenstate of most trotterized operators.
Consequently, achieving an amplitude of $1$ for the fully empty state becomes
difficult.

We further investigate the time evolution of the probabilities in the six-site
system, as depicted in Fig.6. In this plot, we set $\nu=1$ and vary the
trotter step sizes. We plot the probabilities of the fully occupied state and
the empty state, with the solid lines representing the exact results. From the
plot, we observe that while the dynamics effectively capture the departure
from the fully occupied state, achieving a fully empty state is challenging
with the Trotter approximation. However, as we decrease the trotter step size,
the simulation results gradually converge to the exact results.

(a)

(b)

Figure 5: The evolution of the total particle number in the pair annihilation
dynamics for two scenarios: (a) a system with six lattice sites, utilizing a
trotter step size of $\Delta t=1/200$, and (b) a system with seven lattice
sites, employing a trotter step size of $\Delta t=1/50$. In both cases, the
simulation begins with a fully occupied initial condition and the diffusion
constant $D=1$. Figure 6: The evolution of the probabilities of the empty and
fully occupied states in the pair annihilation dynamics with a six-lattice
site system and a fully occupied initial condition: Various trotter step sizes
are tested for the probability of the empty state. The solid lines depict the
exact results, while the scattered markers represent the results obtained from
quantum simulations. The diffusion constant $D=1$ and the reaction rate
$\nu=0.6$.

###  IV.4 Directed Percolation

The last example we explore is a reaction-diffusion system that emulates the
directed percolation phase transition. This system encompasses spontaneous
decay and branching reactions with rates $\nu$ and $\lambda$, respectively.
Additionally, particles are allowed to move freely on the lattice with a rate
$D$. The ‘pseudo’-Hamiltonian governing this system is defined as follows:

$\displaystyle H=-\frac{1}{4}\sum_{\langle
i,j\rangle}\Big{[}2D\sigma_{i}^{x}\sigma_{j}^{x}+2D\sigma_{i}^{y}\sigma_{i}^{y}+2(D+\lambda)\sigma_{i}^{z}\sigma_{j}^{z}$
(24)
$\displaystyle+\lambda\sigma_{i}^{z}\sigma_{j}^{x}+\lambda\sigma_{i}^{x}\sigma_{j}^{z}+i\lambda\sigma_{i}^{z}\sigma_{j}^{y}+i\lambda\sigma_{i}^{y}\sigma_{j}^{z}\Big{]}$
$\displaystyle-\frac{1}{2}\sum_{i}\Big{[}-\nu\sigma_{i}^{z}+(\lambda+\nu)\sigma_{i}^{x}+i(\lambda-\nu)\sigma_{i}^{y}\Big{]}.$

In the regime where $\lambda/\nu\rightarrow 0$, the system converges to an
absorbing stationary state with no particles. Conversely, as $\lambda/\nu\gg
0$, the system enters an active phase where a finite particle density is
sustained. The transition between these phases, occurring for intermediate
values of $\lambda/\nu$, is characterized by an active-absorbing phase
transition. This transition is continuous, and its critical point falls within
the directed percolation universality class.

We simulate a six-site one-dimensional lattice initially populated with two
particles. The results, shown in Fig.7, are obtained with a fixed diffusion
constant $D=1$ and branching rate $\lambda=1$ for simplicity. The scattered
markers represent the results from quantum simulation, while the solid lines
represent the exact results.

As observed, when the decay rate $\nu$ is small, the system maintains a finite
particle density, whereas for large decay rates, the total number of particles
in the system tends towards zero. The results from the quantum simulation
effectively capture the phase transition. However, due to critical slowing-
down near the phase transition, it takes significantly longer to reach the
stationary value, particularly evident for $\nu=0.4$. From the limited data
obtained, we can infer that the critical value $\nu_{c}$ lies between $0.2$
and $0.4$. It’s worth noting that near the critical point, trotter errors tend
to be larger compared to the results obtained deep within the two distinct
phases. This discrepancy arises because deep within the phases, the system
rapidly relaxes to the stationary state. While in the stationary state, the
errors are mostly leakage errors, which are small. However, near the critical
point, the relaxation is slow, leading to the accumulation of trotter errors.
This suggests that a smaller trotter step size is needed when simulating long-
time dynamics, particularly near critical points.

Figure 7: The evolution of particle number in the directed percolation system
with six lattice sites and different decaying rates, $\nu$: The system begins
with an initial condition of two occupied sites. The diffusion constant is set
to $D=1$, and the branching rate is $\lambda=1$. The trotter step size is
chosen as $\Delta t=1/200$. The solid lines depict the exact results, while
the scattered markers represent the results obtained from quantum simulations.

##  V summary and outlook

We explore the digital quantum simulation of classical reaction-diffusion
systems on lattices, focusing on the probability distribution level. To lay
the groundwork, we introduce the spin representation of these systems on
lattices, which provides a convenient framework for simulating them on quantum
computers using qubits. We then introduce the simulation methods of
implementing the dynamics. We start with trotterization, a widely-used
technique for simulating time evolutions that consistently yields robust
results. Next, we introduce the Pauli gadget, which provides a way to
implement time evolution operator of spin systems with universal gate sets.
The Pauli gadget also offers advantages in optimizing quantum circuits. Given
the disparity between the master equation and Schrödinger’s equation,
implementing the dynamics requires handling both unitary and non-unitary
operations. To address this challenge, we employ the probabilistic imaginary
time evolution methods, leveraging Pauli gadgets to reduce the number of
ancilla qubits required. This method aligns well with the nature of spin
representations in the system, particularly facilitating the implementation of
non-unitary operations such as the time evolution operator of the form
$e^{-Ht}$. We also address the pre- and post-processing steps required to
handle the different normalizations between classical probability states and
quantum states. By mapping the normalization constant to a measurement of the
system, we demonstrate that we can circumvent the need for exponential
classical post-processing.

Next, we apply the introduced methods to four distinct examples. Beginning
with the simplest case of a single-site system involving particle generation
and decay, we find that quantum simulation yields highly accurate results even
with a relatively large trotter step size of $\Delta t=1/20$. Moving on, we
investigate free hopping and pair annihilation dynamics on a linear lattice,
focusing on the impact of trotter errors and trotter step size. Our analysis
demonstrates that sufficiently small trotter step sizes provide a reliable
approximation of the exact results, both in terms of probability distributions
and observables.

Finally, we explore a system exhibiting an active-absorbing phase transition.
While quantum simulation accurately captures the phase transition phenomenon,
we observe larger trotter errors near the critical point, necessitating the
use of smaller trotter steps for improved accuracy. These comprehensive
examinations showcase the versatility and efficacy of our quantum simulation
approach across a range of classical reaction-diffusion systems.

One of the primary drawbacks of the method outlined in the paper stems from
the combination of probabilistic imaginary time evolution (PITE) and
Trotterization. As trotter errors accumulate over time, simulating longer
durations necessitates smaller trotter step sizes, leading to an increased
number of trotter steps. However, the probability of measuring the desired
state in the PITE method decays exponentially, demanding exponentially more
realizations of the simulation.

While systems away from critical points typically exhibit rapid relaxation,
critical slowing-down near phase transitions result in significantly longer
relaxation times. In such cases, the combined effects of Trotterization and
PITE pose challenges to the accurate simulation of phase transitions. While
possible methods have been proposed to address the exponential decay of
success probability [94, 95, 90, 96], they still do not fully resolve the
fundamental issue of success probability decay with an increase in trotter
steps. To circumvent this challenge, one possible approach is to replace
Trotterization with the Cartan decomposition [85]. Although requiring
classical preprocessing, the Cartan decomposition offers a way to implement
time evolution operators with fixed circuit depth, thereby avoiding
exponential decay while increasing the number of trotter steps. Alternatively,
variational implementations of non-unitary operators [56, 57, 58, 59, 60] can
be employed to replace PITE. These methods avoid the probabilistic nature of
PITE and instead rely on classical optimizations. Both approaches have the
potential to yield better results for simulating dynamics over longer time
scales, particularly in the study of dynamic phase transitions.

###### Acknowledgements.

The authors express gratitude for valuable discussions with Zhangjie Qin and
proofreading from Uwe C. Täuber. Additionally, the authors acknowledge the
role of AI, particularly ChatGPT, in refining the language of the paper.

## References

  * Janssen [1981] H.-K. Janssen, Zeitschrift für Physik B Condensed Matter 42, 151 (1981).
  * Grassberger [1981] P. Grassberger, in  _Nonlinear Phenomena in Chemical Dynamics: Proceedings of an International Conference, Bordeaux, France, September 7–11, 1981_ (Springer, 1981) p. 262.
  * Cardy and Täuber [1998] J. L. Cardy and U. C. Täuber, Journal of statistical physics 90, 1 (1998).
  * Antal and Droz [2001] T. Antal and M. Droz, Physical Review E 63, 056119 (2001).
  * Deng  _et al._ [2020] S. Deng, W. Li, and U. C. Täuber, Physical Review E 102, 042126 (2020).
  * Howard and Täuber [1997] M. J. Howard and U. C. Täuber, Journal of Physics A: Mathematical and General 30, 7721 (1997).
  * Jensen [1993] I. Jensen, Physical review letters 70, 1465 (1993).
  * Janssen  _et al._ [2004] H.-K. Janssen, F. van Wijland, O. Deloubriere, and U. C. Täuber, Physical Review E 70, 056114 (2004).
  * Nettuno  _et al._ [2024] B. Nettuno, D. Toffenetti, C. Metzl, L. Weigand, F. Raßhofer, R. Swiderski, and E. Frey, arXiv preprint arXiv:2402.06505 (2024).
  * Ben-Naim  _et al._ [2003] E. Ben-Naim, P. L. Krapivsky, and S. Redner, Physica D: nonlinear phenomena 183, 190 (2003).
  * Yao  _et al._ [2023] L. H. Yao, M. Swailem, U. Dobramysl, and U. C. Täuber, Journal of Physics A: Mathematical and Theoretical 56, 225001 (2023).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Reviews of Modern Physics 65, 851 (1993).
  * Butler and Reynolds [2009] T. Butler and D. Reynolds, Physical Review E 79, 032901 (2009).
  * Täuber [2011] U. C. Täuber, Journal of Physics: Conference Series 319, 012019 (2011).
  * Dobramysl  _et al._ [2018] U. Dobramysl, M. Mobilia, M. Pleimling, and U. C. Täuber, Journal of Physics A: Mathematical and Theoretical 51, 063001 (2018).
  * Reichenbach  _et al._ [2007] T. Reichenbach, M. Mobilia, and E. Frey, Nature 448, 1046 (2007).
  * Reichenbach  _et al._ [2008] T. Reichenbach, M. Mobilia, and E. Frey, Journal of Theoretical Biology 254, 368 (2008).
  * Frachebourg  _et al._ [1996] L. Frachebourg, P. L. Krapivsky, and E. Ben-Naim, Physical Review E 54, 6186 (1996).
  * Elton and Nicholson [1942] C. Elton and M. Nicholson, The Journal of Animal Ecology , 215 (1942).
  * Maynard-Smith [1974] J. Maynard-Smith,  _Models in ecology_ (Cambridge university press, 1974).
  * Cheng  _et al._ [2013] L. Y. Cheng, J. D. Murray, and G. M. Odell, in  _Mathematical Topics in Population Biology, Morphogenesis and Neurosciences: Proceedings of an International Symposium held in Kyoto, November 10–15, 1985_ , Vol. 71 (Springer Science & Business Media, 2013) p. 209.
  * Durrett [1999] R. Durrett, SIAM review 41, 677 (1999).
  * Hofbauer and Sigmund [1998] J. Hofbauer and K. Sigmund,  _Evolutionary games and population dynamics_ (Cambridge university press, 1998).
  * Täuber [2012] U. C. Täuber, Journal of Physics A: Mathematical and Theoretical 45, 405002 (2012).
  * Täuber [2014] U. C. Täuber,  _Critical dynamics: a field theory approach to equilibrium and non-equilibrium scaling behavior_ (Cambridge University Press, 2014).
  * Lipowski [1999] A. Lipowski, Physical Review E 60, 5179 (1999).
  * Frey [2010] E. Frey, Physica A: Statistical Mechanics and its Applications 389, 4265 (2010).
  * Henkel and Hinrichsen [2001] M. Henkel and H. Hinrichsen, Journal of Physics A: Mathematical and General 34, 1561 (2001).
  * Henkel  _et al._ [1997] M. Henkel, E. Orlandini, and J. Santos, annals of physics 259, 163 (1997).
  * Alcaraz  _et al._ [1994] F. C. Alcaraz, M. Droz, M. Henkel, and V. Rittenberg, Annals of Physics 230, 250 (1994).
  * Dimentberg [2003] M. F. Dimentberg, in  _IUTAM Symposium on Nonlinear Stochastic Dynamics: Proceedings of the IUTAM Symposium held in Monticello, Illinois, USA, 26–30 August 2002_ (Springer, 2003) p. 307.
  * Swailem and Täuber [2024] M. Swailem and U. C. Täuber, arXiv preprint arXiv:2404.03089 (2024).
  * Carlon  _et al._ [1999] E. Carlon, M. Henkel, and U. Schollwöck, The European Physical Journal B-Condensed Matter and Complex Systems 12, 99 (1999).
  * Helms  _et al._ [2019] P. Helms, U. Ray, and G. K.-L. Chan, Physical Review E 100, 022101 (2019).
  * Merbis  _et al._ [2023] W. Merbis, C. de Mulatier, and P. Corboz, Physical Review E 108, 024303 (2023).
  * Gu and Zhang [2022] J. Gu and F. Zhang, New Journal of Physics 24, 113022 (2022).
  * Kemper  _et al._ [2002] A. Kemper, A. Gendiar, T. Nishino, A. Schadschneider, and J. Zittartz, Journal of Physics A: Mathematical and General 36, 29 (2002).
  * Johnson  _et al._ [2010] T. Johnson, S. Clark, and D. Jaksch, Physical Review E 82, 036702 (2010).
  * Tilly  _et al._ [2022] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, _et al._ , Physics Reports 986, 1 (2022).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, Nature communications 5, 4213 (2014).
  * Tang  _et al._ [2021] H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, arXiv preprint arXiv:1411.4028 (2014).
  * Zhou  _et al._ [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Physical Review X 10, 021067 (2020).
  * Sridhar  _et al._ [2023] V. K. Sridhar, Y. Chen, B. Gard, E. Barnes, and S. E. Economou, arXiv preprint arXiv:2310.09694 (2023).
  * Trabesinger [2012] A. Trabesinger, Nature Physics 8, 263 (2012).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153 (2014).
  * Brown  _et al._ [2010] K. L. Brown, W. J. Munro, and V. M. Kendon, Entropy 12, 2268 (2010).
  * Lloyd [1996] S. Lloyd, Science 273, 1073 (1996).
  * Nielsen  _et al._ [2002] M. A. Nielsen, M. J. Bremner, J. L. Dodd, A. M. Childs, and C. M. Dawson, Physical Review A 66, 022317 (2002).
  * Berry  _et al._ [2015] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Physical review letters 114, 090502 (2015).
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, Physical review letters 118, 010501 (2017).
  * Childs and Wiebe [2012] A. M. Childs and N. Wiebe, arXiv preprint arXiv:1202.5822 (2012).
  * Wang  _et al._ [2011] H. Wang, S. Ashhab, and F. Nori, Physical Review A 83, 062317 (2011).
  * Endo  _et al._ [2020] S. Endo, J. Sun, Y. Li, S. C. Benjamin, and X. Yuan, Physical Review Letters 125, 010501 (2020).
  * Sweke  _et al._ [2015] R. Sweke, I. Sinayskiy, D. Bernard, and F. Petruccione, Physical Review A 91, 062308 (2015).
  * Jones  _et al._ [2019] T. Jones, S. Endo, S. McArdle, X. Yuan, and S. C. Benjamin, Physical Review A 99, 062304 (2019).
  * McArdle  _et al._ [2019] S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, npj Quantum Information 5, 75 (2019).
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. Brandao, and G. K.-L. Chan, Nature Physics 16, 205 (2020).
  * Sun  _et al._ [2021] S.-N. Sun, M. Motta, R. N. Tazhigulov, A. T. Tan, G. K.-L. Chan, and A. J. Minnich, PRX Quantum 2, 010317 (2021).
  * Nishi  _et al._ [2021] H. Nishi, T. Kosugi, and Y.-i. Matsushita, npj Quantum Information 7, 85 (2021).
  * Leadbeater  _et al._ [2023] C. Leadbeater, N. Fitzpatrick, D. M. Ramo, and A. J. Thom, arXiv preprint arXiv:2304.07917 (2023).
  * Cao  _et al._ [2013] Y. Cao, A. Papageorgiou, I. Petras, J. Traub, and S. Kais, New Journal of Physics 15, 013021 (2013).
  * Berry [2014] D. W. Berry, Journal of Physics A: Mathematical and Theoretical 47, 105301 (2014).
  * Montanaro and Pallister [2016] A. Montanaro and S. Pallister, Physical Review A 93, 032324 (2016).
  * Linden  _et al._ [2022] N. Linden, A. Montanaro, and C. Shao, Communications in Mathematical Physics 395, 601 (2022).
  * Kubo  _et al._ [2021] K. Kubo, Y. O. Nakagawa, S. Endo, and S. Nagayama, Physical Review A 103, 052425 (2021).
  * Banchi [2023] L. Banchi, Machine Learning: Science and Technology (2023).
  * Ghafari  _et al._ [2019] F. Ghafari, N. Tischler, J. Thompson, M. Gu, L. K. Shalm, V. B. Verma, S. W. Nam, R. B. Patel, H. M. Wiseman, and G. J. Pryde, Physical Review X 9, 041013 (2019).
  * Korzekwa and Lostaglio [2021] K. Korzekwa and M. Lostaglio, Physical Review X 11, 021019 (2021).
  * Kiss  _et al._ [2023] O. Kiss, M. Grossi, and A. Roggero, Quantum 7, 977 (2023).
  * Blank  _et al._ [2021] C. Blank, D. K. Park, and F. Petruccione, npj Quantum Information 7, 126 (2021).
  * Kalantzis  _et al._ [2024] V. Kalantzis, M. S. Squillante, and S. Ubaru, On quantum algorithms for efficient solutions of general classes of structured markov processes (2024), arXiv:2404.17959 [quant-ph] .
  * Yepez [1998] J. Yepez, in  _NASA International Conference on Quantum Computing and Quantum Communications_ (Springer, 1998) p. 34.
  * Gaitan [2020] F. Gaitan, npj Quantum Information 6, 61 (2020).
  * Steijl [2019] R. Steijl, Advances in quantum communication and information 31 (2019).
  * Kapil  _et al._ [2018] M. Kapil, B. K. Behera, and P. K. Panigrahi, arXiv preprint arXiv:1807.00521 (2018).
  * Jafferis  _et al._ [2022] D. Jafferis, A. Zlokapa, J. D. Lykken, D. K. Kolchmeyer, S. I. Davis, N. Lauk, H. Neven, and M. Spiropulu, Nature 612, 51 (2022).
  * Kamakari  _et al._ [2022] H. Kamakari, S.-N. Sun, M. Motta, and A. J. Minnich, PRX Quantum 3, 010320 (2022).
  * Pagano  _et al._ [2020] G. Pagano, A. Bapat, P. Becker, K. S. Collins, A. De, P. W. Hess, H. B. Kaplan, A. Kyprianidis, W. L. Tan, C. Baldwin, _et al._ , Proceedings of the National Academy of Sciences 117, 25396 (2020).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, _et al._ , Science 374, 1242 (2021).
  * Berry  _et al._ [2014] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, in  _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_ (2014) p. 283.
  * Childs  _et al._ [2018] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences 115, 9456 (2018).
  * Wen [2004] X.-G. Wen,  _Quantum field theory of many-body systems: From the origin of sound to an origin of light and electrons_ (OUP Oxford, 2004).
  * Childs  _et al._ [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Physical Review X 11, 011020 (2021).
  * Kökcü  _et al._ [2022] E. Kökcü, T. Steckmann, Y. Wang, J. Freericks, E. F. Dumitrescu, and A. F. Kemper, Physical Review Letters 129, 070501 (2022).
  * Van Den Berg and Temme [2020] E. Van Den Berg and K. Temme, Quantum 4, 322 (2020).
  * Mukhopadhyay  _et al._ [2023] P. Mukhopadhyay, N. Wiebe, and H. T. Zhang, npj Quantum Information 9, 31 (2023).
  * Gingrich and Williams [2004] R. M. Gingrich and C. P. Williams, in  _ACM International Conference Proceeding Series_ , Vol. 58 (Citeseer, 2004) p. 1.
  * Terashima and Ueda [2005] H. Terashima and M. Ueda, International Journal of Quantum Information 3, 633 (2005).
  * Nishi  _et al._ [2022] H. Nishi, T. Kosugi, Y. Nishiya, and Y. ichiro Matsushita, Acceleration of probabilistic imaginary-time evolution method combined with quantum amplitude amplification (2022), arXiv:2212.13816 [quant-ph] .
  * Liu  _et al._ [2021] T. Liu, J.-G. Liu, and H. Fan, Quantum Information Processing 20, 204 (2021).
  * Lin  _et al._ [2021] S.-H. Lin, R. Dilip, A. G. Green, A. Smith, and F. Pollmann, PRX Quantum 2, 010342 (2021).
  * Silva  _et al._ [2023] T. L. Silva, M. M. Taddei, S. Carrazza, and L. Aolita, Scientific Reports 13, 18258 (2023).
  * Grover [1996] L. K. Grover, in  _Proceedings of the twenty-eighth annual ACM symposium on Theory of computing_ (1996) p. 212.
  * Grover [2005] L. K. Grover, Physical Review Letters 95, 150501 (2005).
  * Chen  _et al._ [2024] Y. Chen, I. Fulginiti, and C. B. Mendl, Reducing mid-circuit measurements via probabilistic circuits (2024), arXiv:2405.13747 [quant-ph] .