# Measurement of the electron-hole pair creation energy in a 4H-SiC p-n diode

Andreas Gsponer<EMAIL_ADDRESS>Matthias Knopf  Philipp Gaggl
Jürgen Burin  Simon Waid  Thomas Bergauer  Institute of High Energy Physics,
Austrian Academy of Sciences, Vienna, Austria

###### Abstract

For 4H silicon carbide (4H-SiC), the values for the electron-hole pair
creation energy $\epsilon_{\text{i}}$ published in literature vary
significantly. This work presents an experimental determination of
$\epsilon_{\text{i}}$ for a $50$ $\mu$m 4H-SiC p-n diode designed for particle
detection in high energy physics. The detector response was measured for
$\alpha$ particles between $4.2$ MeV and $5.6$ MeV for 4H-SiC and a silicon
reference device. Different $\alpha$ energies were obtained by using multiple
nuclides and varying the effective air gap between the $\alpha$ source and the
detector. The energy deposited in the detectors was determined using a Monte
Carlo simulation taking into account the device cross-section. A linear fit of
the detector response to the deposited energy yields
$\epsilon_{\text{i}}=(7.83\pm 0.02)\;\text{eV}$, which agrees well with the
most recent literature.

###### keywords:

Silicon carbide , 4H-SiC , Wide-bandgap detector , Electron-pair creation
energy , planar diode

††journal: NIM A

##  1 Introduction

4H-Silicon carbide (4H-SiC) is gaining popularity as a material for detector
development. It offers multiple advantages over the ubiquitous silicon (Si)
detectors currently in use. The wide bandgap of 4H-SiC and large displacement
threshold [1] allows for high-temperature operation [2] and enables
potentially higher resilience against radiation damage (especially relevant
for future high luminosity experiments). The high breakdown voltage and
carrier saturation velocity of 4H-SiC promises to enable enhanced timing
performances compared to Si [3, 4].

For the design of 4H-SiC detectors and their readout electronics, a proper
understanding of 4H-SiC material properties is essential. Among these are the
ionization energy $\epsilon_{\text{i}}$ (energy needed to create an electron-
hole pair) and the intrinsic energy resolution described by the Fano factor F.
In the literature, there is still a considerable disagreement over their
respective values, likely also due to the ongoing considerable improvements in
material quality [1].

In this paper, we thus provide an overview of 4H-SiC values for
$\epsilon_{\text{i}}$ and $F$ available in the literature. We further present
measurements on 4H-SiC- and a Si-based particle detectors, both designed
primarily for high-energy physics (HEP) applications, to determine
$\epsilon_{\text{i}}$ of a 4H-SiC p-n diode. Differences in the detector
cross-section are accounted for using a precise GATE [5] simulation. By
controlling the air pressure during the measurement (from vacuum to ambient
air pressure), a range of different energies and distributions can be accessed
using the same $\alpha$ source. For extracting $\epsilon_{\text{i}}$, we
finally compare the signals for 4H-SiC and Si and scale the ratio with the
well-known ionization energy of the latter.

###  1.1 Ionization Energy and Fano Factor

Table 1: Previously published values of $\epsilon_{\text{i,SiC}}$ for 4H-SiC. For measurements that compare to silicon, the value of $\epsilon_{\text{i,Si}}$ is indicated. The reference labeled “*” refers to the results obtained in this work. Ref. | Year |  $\epsilon_{\text{i,SiC}}$ [eV] | Radiation | Method | Device Type |  $\epsilon_{\text{i, Si}}$ [eV]
---|---|---|---|---|---|---
* | 2023 | $7.83\pm 0.02$ | $\alpha$ | Comp. to Si | Diode | 3.62
[6] | 2013 | $7.82\pm 0.02$ | $\alpha$ | Charge inj. | Schottky | -
[7] | 2013 | 7.28 | $\alpha$ | Charge inj. | Schottky | -
[8] | 2006 | 5.05 | $e^{-}$ | DC Current | Diode | -
[9] | 2005 | $7.78\pm 0.05$ | $\alpha$ | Comp. to Si | Schottky | 3.62
[9] | 2005 | $7.79\pm 0.09$ | protons | Comp. to Si | Schottky | 3.64
[10] | 2005 | 7.6 |  X-Rays (^241Am) | Comp. to Si | Diode | 3.60
[11] | 2004 | 8.6 | $\alpha$ | Comp. to Si | Diode | 3.62
[12] | 2004 | 7.71 | $\alpha$ | Comp. to Si | Schottky | 3.62
[13] | 2003 | 7.8 | X-rays (SEM) | Comp. to Si | Diode | 3.67
Table 2: Values of $F_{\text{SiC}}$ for 4H-SiC used in published literature. Ref. | Year | $F$ | Radiation | Method | Device Type
---|---|---|---|---|---
[7] | 2013 | 0.128 | $\alpha$ |  Klein model [14] | Schottky
[15] | 2011 | 0.1 | X-rays | Calculated from linewidth | Schottky
[10] | 2005 | $\leq 0.04$ | X-rays | Upper limit from linewidth | Diode
[13] | 2003 | 0.12 | - | Hypothesized | -

The average energy a single primary particle has to expend to create one
electron-hole pair on its passage through matter is referred to as the
_ionization energy_ $\epsilon_{\text{i}}$. It can be used to convert the
incident energy of a primary particle to an amount of generated charge carrier
pairs $n_{\text{e-h}}$, provided that the particle stops in the medium, thus
transferring its total kinetic energy In an ideal detector, all events along a
particle track are considered independent, and the total number of electron-
hole pairs is equal to $n_{\text{e-h}}=E/\epsilon_{\text{i}}$, where $E$ is
the total deposited energy. The stopping of particles in matter is subject to
energy straggling and exhibits statistical fluctuations in the number of
individual processes and energy losses. This process may be described by the
Poisson statistic and leads to signal fluctuations on the order of
$\sqrt{n_{\text{e-h}}}$. However, due to the fact that the individual
particle-particle interactions are not statistically independent [16], the
energy resolution is improved by a factor $\sqrt{F}$. This _Fano factor_ $F$
is a measure of the intrinsic energy resolution for a given material due to
fluctuations in the primary signal formation process. It is defined as the
ratio of the observed variance $\sigma^{2}$ to the variance of the Poisson
statistic $\sqrt{n_{\text{e-h}}}$

$F=\frac{\text{observed variance}}{\text{Poisson
variance}}=\frac{\sigma^{2}}{E/\epsilon_{\text{i}}}.$ (1)

With this, the standard deviation in terms of energy
$\sigma_{\text{Fano}}=\sigma/\epsilon_{\text{i}}$ can be written in the usual
form:

$\sigma_{\text{Fano}}=\sqrt{\epsilon_{\text{i}}FE}$ (2)

In practice, however, there are also other effects contributing to the
variation of the measured energy resolution $\sigma_{\text{meas.}}$, i.e.,

$\sigma_{\text{meas.}}^{2}=\sigma_{\text{Fano}}^{2}+\sigma_{\text{noise}}^{2}+\sigma_{\text{absorbed}}^{2}$
(3)

The contribution $\sigma_{\text{noise}}$ corresponds to the noise of the
readout electronics and $\sigma_{\text{absorbed}}$ to energy straggling caused
by material that covers the sensitive volume of the detector [7]. Losses due
to the variation in charge collection can be minimized by using high-quality
materials and adequate bias voltages.

In the literature, a wide range of ionization energies are reported for
4H-SiC, with values ranging between $5.05\text{\,}\mathrm{eV}$ [8] and
$8.6\text{\,}\mathrm{eV}$ [11]. A summary is shown in Table 1. However, more
recent measurements have converged in the range of
$7.6\text{\,}\mathrm{eV}7.8\text{\,}\mathrm{eV}$ [7, 12, 10, 13]. Different
methods have been used, but most frequently, $\epsilon_{\text{i}}$ was
determined via a comparison to a silicon detector. This method is described in
more detail in the next section, 1.2. As the ionization energies obtained by
this method depend directly on the $\epsilon_{\text{i}}$ of silicon,
$\epsilon_{\text{i,Si}}$, the values used in each publication are listed in
the table as well. Due to the difficulty of measuring the Fano factor, only a
few published values are available for 4H-SiC. An overview is given in table
2. Estimations have been made from the Fano factor of silicon [13], and using
Klein’s model [14, 7]. Experimental data is available from X-ray measurements,
with one publication observing no significant statistical broadening of the
energy resolution (and estimating $F_{\text{SiC}}<$0.04$$ [10]), and another
finding $F_{\text{SiC}}=$0.100$$ [15].

###  1.2 Comparison Method with Silicon

To determine $\epsilon_{\text{i}}$ based on a comparison to Si, a
spectroscopic measurement is performed in the same experimental setup for a Si
and a 4H-SiC detector. In most cases, $\alpha$ sources are used for this
purpose, as $\alpha$ particles are quickly absorbed even in thin detectors
(with a range of about $26\text{\,}\mathrm{\SIUnitSymbolMicro m}$ for ^241Am
in Si [17]). In an ideal experimental setup, the ionization energy can then be
calculated from the ratio of the obtained signals ($\mu_{\text{Si}}$ and
$\mu_{\text{SiC}}$), which is directly proportional to the number of generated
charge carriers $n_{\text{e-h}}$:

$\frac{\mu_{\text{Si}}}{\mu_{\text{SiC}}}=\frac{E_{\alpha}\epsilon_{\text{i,SiC}}}{E_{\alpha}\epsilon_{\text{i,Si}}}\Rightarrow\epsilon_{\text{i,SiC}}=\frac{\mu_{\text{Si}}}{\mu_{\text{SiC}}}\cdot\epsilon_{\text{i,Si}}.$
(4)

Here, $\epsilon_{\text{i,Si}}$ and $\epsilon_{\text{i,SiC}}$ denote the
ionization energies of Si and 4H-SiC respectively, and $E_{\alpha}$ is the
energy deposited by the $\alpha$ particles. The main advantage of this method
is that no absolute charge calibration of the electronics is necessary [7].

However, there are a couple of caveats: This method assumes a charge
collection efficiency (CCE) of $100\text{\,}\mathrm{\char 37\relax}$ and an
identical amplifier response for both detectors. The former can be achieved by
operating the detector with voltages above the full depletion voltage, while
the latter can be minimized by using devices with similar capacitance or
introducing a correction factor [9]. A major source of uncertainties is the
energy loss in the detector layers above the sensitive volume (passivation and
metal layers). If there is a difference in these layers among devices, this
needs to be considered as well, using stopping power calculations [9] or Monte
Carlo simulations [7].

##  2 Materials and Methods

The employed 4H-SiC detector was a planar p-n diode developed and manufactured
by IMB-CNM-CSIC [18]. This detector type had previously been investigated
using UV-TCT [19], $\alpha$ particles [20] and proton beams [21, 22], before
and after irradiation with neutrons.

Figure 1: Si (left) and 4H-SiC (right) detectors mounted on ceramic PCBs.

This detector featured a $50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ epitaxial
layer on a 4H-SiC substrate, with an area of $3\times 3$ mm and a resistivity
of $20\text{\,}\mathrm{\SIUnitSymbolOhm}\text{\,}\mathrm{cm}$. Full depletion
was attained at $300\text{\,}\mathrm{V}$. However, for the $\alpha$ source a
charge collection efficiency of $100\text{\,}\mathrm{\char 37\relax}$ was
already reached at $100\text{\,}\mathrm{V}$, due to the low penetration depth.
The active 4H-SiC volume was covered by a metal layer
($1020\text{\,}\mathrm{nm}$) consisting of titanium, aluminum, and nickel as
well as a passivation layer (SiO_2, Si_3N_4) [18, 23]. The thickness of the
layers has been measured using transmission electron microscopy (TEM) and
corresponds to an equivalent aluminum thickness of
$1.72\text{\,}\mathrm{\SIUnitSymbolMicro m}$. A picture of the 4H-SiC detector
and the reference Si device is shown in Figure 1.

The silicon device was a planar n-type Si diode from the CMS production [24].
This diode had an active thickness of around
$300\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and a measured full depletion
voltage of $180\text{\,}\mathrm{V}$. The metal layer consisted of only
aluminum (albeit thicker than for the 4H-SiC sample), with a passivation layer
very similar to the 4H-SiC device, for an equivalent aluminum thickness of
$2.23\text{\,}\mathrm{\SIUnitSymbolMicro m}$. During all measurements, the Si
and 4H-SiC detectors were operated using a reverse bias of
$400\text{\,}\mathrm{V}$ to ensure full depletion and thus a
$100\text{\,}\mathrm{\char 37\relax}$ charge collection efficiency. The alpha
spectra were measured using an Eckert & Ziegler QCRB2500 spectrometric mixed
nuclide source. The radionuclides present in the source were ^239Pu, ^241Am
and ^244Cm, with an activity of $1\text{\,}\mathrm{k}\mathrm{B}\mathrm{q}$ per
isotope. The main decay energies of the employed sources are almost
equidistant, at $5.157\text{\,}\mathrm{MeV}$, $5.483\text{\,}\mathrm{MeV}$ and
$5.805\text{\,}\mathrm{MeV}$ [25]. The radionuclides were concentrated on a
$7\text{\,}\mathrm{m}\mathrm{m}$ disk with a thickness of about
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$, located behind a steel collimator.

The detectors were mounted and wire-bonded on passive ceramic PCBs (Figure 1).
A custom-made holder was designed, which allows for a constant source-detector
distance of $6\text{\,}\mathrm{m}\mathrm{m}$. The entire setup was wrapped
with aluminum foil and situated inside a stainless steel vacuum chamber to
ensure RF shielding. The spectroscopic measurements were performed with a
readout chain comprising of a Cividec Cx-L shaping charge sensitive amplifier
(CSA). The pulse heights were measured and histogrammed using a Rhode &
Schwarz RTO6 oscilloscope in HD mode ($16\text{\,}\mathrm{bit}$ resolution at
an analog bandwidth of $100\text{\,}\mathrm{MHz}$). A Keithley 2470 source
meter was used to apply the reverse bias via a bias-T integrated into the CSA.

When alpha particles travel through the air, they scatter and quickly lose
energy (with an average range in the order of $\mathrm{cm}$). Small changes in
the source-detector distance can lead to significant uncertainties if the
experimental geometry is not reproducible. In order to eliminate this effect,
measurements can be performed in vacuum. In this work, however, we performed
measurements in vacuum and at different air pressures to obtain results for
multiple $\alpha$-particle energies. This additionally allows us to verify the
simulation results of the energy lost in the passivation and metal layers. In
order to perform measurements at different air pressures, a Pfeiffer EVR116
valve and a Pfeiffer RPT200 pressure gauge were used to actively control the
pressure inside the setup using a PID loop. All measurements were performed at
room temperature ($25\text{\,}\mathrm{\SIUnitSymbolDegree}$).

As the individual decay energies of the nuclides could not be separated in the
data (due to the aforementioned straggling), one single fit was used per
nuclide. The individual decay lines (and their probability) introduce an
asymmetric energy distribution, so a skew normal distribution [26] was used
for the fits. The probability density function $\phi(x)$ of this distribution
is given by

$\phi(x)=\frac{2}{\omega\sqrt{2\pi}}e^{-\frac{\left(x-\xi\right)^{2}}{2\omega^{2}}}\int_{-\infty}^{\alpha\left(\frac{x-\xi}{\omega}\right)}\frac{1}{\sqrt{2\pi}}e^{-\frac{t^{2}}{2}}\mathrm{d}t,$
(5)

with $\xi$ being the location, $\omega$ the scale and $\alpha$ the shape
parameters respectively. The mean $\mu$ and variance $\sigma^{2}$ of this
distribution are given by

$\mu=\xi+\omega\frac{\alpha}{\sqrt{1+\alpha^{2}}}\sqrt{\frac{2}{\pi}}$ (6)

and

$\sigma^{2}=\omega^{2}\left(1-\frac{2\alpha^{2}}{\left(1+\alpha^{2}\right)\pi}\right).$
(7)

###  2.1 Simulation Results

A GATE [5] simulation was used to determine the energy spectrum expected in
the detectors (Si and 4H-SiC) as a function of the air pressure in the vacuum
chamber. For each nuclide and decay line, 5 million events were simulated, and
the total energy deposited inside the detector was histogrammed. Figure 2
shows the average energy per nuclide for both detector types. As expected, the
energy lost in the air gap increases with the air pressure. An energy range of
around $0.5\text{\,}\mathrm{MeV}$ can be accessed by tuning the pressure. Due
to the thicker metal layer of the Si device, alpha-particles of about
$150\text{\,}\mathrm{keV}$ less energy reaches the sensitive volume than for
the 4H-SiC device.

Figure 2: Energy deposited in detector $\text{E}_{\alpha}$ as a function of
air pressure and nuclide for Si (black) and 4H-SiC (red).

##  3 Results

Measurements were performed in air at pressures between
$0.3\text{\,}\mathrm{mbar}$ and $1000\text{\,}\mathrm{mbar}$. For each
pressure setting, 25k events were recorded for the Si and 4H-SiC detectors
individually. Figure 3 shows a comparison of the measured and simulated
spectra at $0.3\text{\,}\mathrm{mbar}$.

(a) SiC spectrum

(b) Si spectrum

Figure 3: Histograms of the measured alpha spectrum at a vacuum of
$0.3\text{\,}\mathrm{mbar}$ for the 4H-SiC (red) and Si (grey) detectors. The
simulated energy deposition is indicated by a solid line, and the measured
data is shown as dots.

Compared to the expected spectrum, the line width of the experimental data is
significantly broadened. This can mainly be attributed due to a non-ideal RF
shielding of the experimental setup. The noise of the CSA was measured to be
$6.2\text{\,}\mathrm{k\text{$e$}}$, which corresponds to a 1-$\sigma$
broadening of around $22\text{\,}\mathrm{keV}$ for the Si and
$50\text{\,}\mathrm{keV}$ for the 4H-SiC detectors. As this broadening is much
larger than the Fano contribution, estimated to below a standard deviation of
$2\text{\,}\mathrm{keV}$ (compare with [7]), these spectra do not allow for an
estimation of the Fano factor for 4H-SiC.

Figure 4 shows the fitted center positions plotted versus the energy deposited
in the detector (depicted in Figure 2). Linear fits were performed for the Si
and 4H-SiC data, resulting in a slope of $a_{\text{Si}}=$0.398\pm
0.001\text{\,}\mathrm{V}\text{\,}{\mathrm{MeV}}^{-1}$$ and
$a_{\text{SiC}}=$0.184\pm
0.002\text{\,}\mathrm{V}\text{\,}{\mathrm{MeV}}^{-1}$$.

Figure 4: Detector signal versus energy in the detector. The color of the
measured data points corresponds to the pressure inside of the vacuum chamber.

Taking the ratio of the slopes (compare with equation 4) and using a silicon
ionization energy $\epsilon_{\text{i,Si}}=$$3.62\text{\,}\mathrm{eV}$ [27],
the ionization energy for 4H-SiC is obtained:

$\epsilon_{\text{i,SiC}}=$7.83\pm 0.02\text{\,}\mathrm{eV}$$

##  4 Discussion

The ionization energy of 4H silicon carbide has been determined using a p-n
diode designed for high energy physics applications. The spectroscopic
measurement leveraged an $\alpha$ source and an air gap between the source and
detector at different air pressures, resulting in varying energy deposition in
the detector. By comparing the spectrum obtained for the 4H-SiC device with a
silicon reference, an ioniziation energy of $\epsilon_{\text{i,SiC}}=$7.83\pm
0.02\text{\,}\mathrm{eV}$$ was obtained. This value compares well with the
most recent results in the literature (see summary table 1). The energy
spectra measured for both the Si and 4H-SiC devices showed a larger linewidth
than expected from Monte-Carlo simulations. This effect can mainly be
attributed to non-ideal RF shielding and large area diodes (and the associated
capacitance), which increased the noise equivalent of the amplifier to
$6.2\text{\,}\mathrm{keV}$. Due to these large linewidths, no estimation of
the Fano factor for 4H-SiC was possible.

Future measurements will aim to reduce the electronic noise in order to
determine the effect of Fano noise. Additionally, devices with thinner
passivation and metal layers promise to reduce the energy straggling in the
entrance window. However, if these layers are precisely known, their effect
could also be compensated by Monte-Carlo simulations. Finally, an absolute
calibration using a precision pulser could be performed for the amplifier,
similar to [7], which can serve as an additional verification of the results
obtained in this paper.

## Acknowledgments

This project has received funding from the Austrian Research Promotion Agency
FFG, grant number 883652. Production and development of the 4H-SiC samples was
supported by the Spanish State Research Agency (AEI) and the European Regional
Development Fund (ERDF), ref. RTC-2017-6369-3.

## References

  * [1] M. De Napoli, Sic detectors: A review on the use of silicon carbide as radiation detection material, Frontiers in Physics 10 (2022).  doi:10.3389/fphy.2022.898833.
  * [2] F. Nava, G. Bertuccio, A. Cavallini, E. Vittone, Silicon carbide and its use as a radiation detector material, Measurement Science and Technology 19 (10) (2008) 102001.  doi:10.1088/0957-0233/19/10/102001.
  * [3] T. Yang, Y. Tan, Q. Liu, S. Xiao, K. Liu, J. Zhang, R. Kiuchi, M. Zhao, X. Zhang, C. Wang, B. Wu, J. Lin, W. Song, H. Lu, X. Shi, Time resolution of the 4h-SiC PIN detector, Frontiers in Physics 10 (mar 2022).  doi:10.3389/fphy.2022.718071.
  * [4] P. Barletta, M. Cerullo, C. Haber, S. E. Holland, J. Muth, B. Sekely, Fast timing with silicon carbide low gain avalanche detectors, type: article (2022).  arXiv:2203.08554[physics], doi:10.48550/arXiv.2203.08554.
  * [5] S. Jan, G. Santin, D. Strul, S. Staelens, K. Assié, D. Autret, S. Avner, R. Barbier, M. Bardiès, P. M. Bloomfield, D. Brasse, V. Breton, P. Bruyndonckx, I. Buvat, A. F. Chatziioannou, Y. Choi, Y. H. Chung, C. Comtat, D. Donnarieix, L. Ferrer, S. J. Glick, C. J. Groiselle, D. Guez, P.-F. Honore, S. Kerhoas-Cavata, A. S. Kirov, V. Kohli, M. Koole, M. Krieguer, D. J. van der Laan, F. Lamare, G. Largeron, C. Lartizien, D. Lazaro, M. C. Maas, L. Maigne, F. Mayet, F. Melot, C. Merheb, E. Pennacchio, J. Perez, U. Pietrzyk, F. R. Rannou, M. Rey, D. R. Schaart, C. R. Schmidtlein, L. Simon, T. Y. Song, J.-M. Vieira, D. Visvikis, R. V. de Walle, E. Wieërs, C. Morel, GATE: a simulation toolkit for PET and SPECT, Physics in Medicine and Biology 49 (19) (2004) 4543–4561.  doi:10.1088/0031-9155/49/19/007.
  * [6] T. R. Garcia, A. Kumar, B. Reinke, T. E. Blue, W. Windl, Electron-hole pair generation in SiC high-temperature alpha particle detectors, Applied Physics Letters 103 (15) (2013) 152108.  doi:10.1063/1.4824774.
  * [7] S. K. Chaudhuri, K. J. Zavalla, K. C. Mandal, Experimental determination of electron-hole pair creation energy in 4H-SiC epitaxial layer: An absolute calibration approach, Applied Physics Letters 102 (3) (2013) 031109.  doi:10.1063/1.4776703.
  * [8] M. V. S. Chandrashekhar, C. I. Thomas, M. G. Spencer, Measurement of the mean electron-hole pair ionization energy in 4h sic, Applied Physics Letters 89 (4) (2006) 042113.  doi:10.1063/1.2243799.
  * [9] A. Lo Giudice, F. Fizzotti, C. Manfredotti, E. Vittone, F. Nava, Average energy dissipated by mega-electron-volt hydrogen and helium ions per electron-hole pair generation in 4h-sic, Applied Physics Letters 87 (22) (2005) 222105.  doi:10.1063/1.2135507.
  * [10] B. Phlips, K. Hobart, F. Kub, R. Stahlbush, M. Das, B. Hull, G. De Geronimo, P. O’Connor, Silicon carbide pin diodes as radiation detectors, in: IEEE Nuclear Science Symposium Conference Record, 2005, Vol. 3, 2005, pp. 1236–1239.  doi:10.1109/NSSMIC.2005.1596542.
  * [11] A. Lebedev, A. Ivanov, N. Strokan, Radiation resistance of SiC and nuclear-radiation detectors based on SiC films, Semiconductors 38 (2) (2004) 125–147.
  * [12] A. Ivanov, E. Kalinina, A. Konstantininov, N. Onushkin, G.A.and Strokan, G. Kholuyanov, A. Hallén, High-resolution short range ion detectors based on 4h-sic films, Technical Physics Letters 30 (2004) 575–577.  doi:10.1134/1.1783406.
  * [13] G. Bertuccio, R. Casiraghi, Study of silicon carbide for x-ray detection and spectroscopy, IEEE Transactions on Nuclear Science 50 (1) (2003) 175–185.  doi:10.1109/TNS.2003.807855.
  * [14] C. A. Klein, Bandgap dependence and related features of radiation ionization energies in semiconductors, Journal of Applied Physics 39 (4) (1968) 2029–2038.  doi:10.1063/1.1656484.
  * [15] G. Bertuccio, S. Caccia, D. Puglisi, D. Macera, Advances in silicon carbide x-ray detectors, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 652 (1) (2011) 193–196, symposium on Radiation Measurements and Applications (SORMA) XII 2010.  doi:10.1016/j.nima.2010.08.046.
  * [16] U. Fano, Ionization yield of radiations. ii. the fluctuations of the number of ions, Phys. Rev. 72 (1947) 26–29.  doi:10.1103/PhysRev.72.26.
URL https://link.aps.org/doi/10.1103/PhysRev.72.26

  * [17] M. J. Berger, J. Coursey, M. Zucker, J. Chang, Stopping-power & range tables for electrons, protons, and helium ions, NIST Standard Reference Database 124 (2017).
  * [18] J. M. Rafí, G. Pellegrini, P. Godignon, S. O. Ugobono, G. Rius, I. Tsunoda, M. Yoneoka, K. Takakura, G. Kramberger, M. Moll, Electron, neutron, and proton irradiation effects on sic radiation detectors, IEEE Transactions on Nuclear Science 67 (12) (2020) 2481–2489.  doi:10.1109/TNS.2020.3029730.
  * [19] P. Gaggl, T. Bergauer, M. Göbel, R. Thalmeier, M. Villa, S. Waid, Charge collection efficiency study on neutron-irradiated planar silicon carbide diodes via uv-tct, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 1040 (2022) 167218.  doi:10.1016/j.nima.2022.167218.
  * [20] P. Gaggl, A. Gsponer, R. Thalmeier, S. Waid, G. Pellegrini, P. Godignon, J. Rafí, T. Bergauer, Performance of neutron-irradiated 4h-silicon carbide diodes subjected to alpha radiation, Journal of Instrumentation 18 (01) (2023) C01042.  doi:10.1088/1748-0221/18/01/C01042.
  * [21] M. Christanell, M. Tomaschek, T. Bergauer, 4h-silicon carbide as particle detector for high-intensity ion beams, Journal of Instrumentation 17 (01) (2022) C01060.  doi:10.1088/1748-0221/17/01/c01060.
  * [22] A. Gsponer, P. Gaggl, J. Maier, R. Thalmeier, S. E. Waid, T. Bergauer, Neutron radiation induced effects in 4h-sic pin diodes (2023).  doi:10.48550/ARXIV.2310.02047.
  * [23] J. Rafí, G. Pellegrini, P. Godignon, D. Quirion, S. Hidalgo, O. Matilla, A. Fontserè, B. Molas, K. Takakura, I. Tsunoda, M. Yoneoka, D. Pothin, P. Fajardo, Four-quadrant silicon and silicon carbide photodiodes for beam position monitor applications: electrical characterization and electron irradiation effects, Journal of Instrumentation 13 (01) (2018) C01045.  doi:10.1088/1748-0221/13/01/C01045.
  * [24] J.-L. Agram, M. Angarano, S. Assouak, T. Bergauer, G. Bilei, L. Borrello, M. Brianzi, C. Civinini, A. Dierlamm, N. Dinu, N. Demaria, L. Feld, E. Focardi, J.-C. Fontaine, E. Forton, A. Furgeri, G. Gregoire, F. Hartmann, A. Honma, P. Juillot, D. Kartashov, M. Krammer, A. Macchiolo, M. Mannelli, A. Messineo, E. Migliore, O. Militaru, C. Piasecki, R. Santinelli, D. Sentenac, L. Servoli, A. Starodumov, G. Tonelli, J. Wang, The silicon sensors for the compact muon solenoid tracker—design and qualification procedure, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 517 (1-3) (2004) 77–93.  doi:10.1016/j.nima.2003.08.175.
  * [25] A. Sonzogni, Nudat 3, National Nuclear Data Center (NNDC) - Brookhaven National Laboratory (2023).
  * [26] A. Azzalini, A. Capitanio, The Skew-Normal and Related Families, Cambridge University Press, 2013.  doi:10.1017/cbo9781139248891.
  * [27] M. Mazziotta, Electron–hole pair creation energy and fano factor temperature dependence in silicon, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 584 (2-3) (2008) 436–439.  doi:10.1016/j.nima.2007.10.043.


# Quantum Finite Volume Method for Computational Fluid Dynamics
with Classical Input and Output

Zhao-Yun Chen  Cheng Xue  Si-Ming Chen  Bing-Han Lu  Yu-Chun Wu
<EMAIL_ADDRESS>Key Laboratory of Quantum Information, CAS  University of
Science and Technology of China  Ju-Chun Ding  Sheng-Hong Huang
<EMAIL_ADDRESS>Department of Modern Mechanics, USTC  Guo-Ping Guo  Key
Laboratory of Quantum Information, CAS  University of Science and Technology
of China  Origin Quantum Computing, Hefei<EMAIL_ADDRESS>
###### Abstract

Computational fluid dynamics (CFD) is a branch of fluid mechanics that uses
numerical methods to solve fluid flows. The finite volume method (FVM) is an
important one. In FVM, space is discretized to many grid cells. When the
number of grid cells grows, massive computing resources are needed
correspondingly. Recently, quantum computing has been proven to outperform a
classical computer on specific computational tasks. However, the quantum CFD
(QCFD) solver remains a challenge because the conversion between the classical
and quantum data would become the bottleneck for the time complexity. Here we
propose a QCFD solver with exponential speedup over classical counterparts and
focus on how a quantum computer handles classical input and output. By
utilizing quantum random access memory, the algorithm realizes sublinear time
at every iteration step. The QCFD solver could allow new frontiers in the CFD
area by allowing a finer mesh and faster calculation.

##  I Introduction

Computational fluid dynamics (CFD) is the area that utilizes numerical methods
to obtain the physical properties of fluids. It has many applications, such as
aid in designing aircraft or automobile. CFD is often related to solving a
series of partial differential equations (PDEs) and can compute the evolution
of physical characteristics of fluid at a given space, including density,
momentum, and energy. These characteristics would provide us essential
references for the properties of the fluid in the computational space. There
are three typical physical governing equations of the CFD: Navier–Stokes (NS),
Euler, and Reynold–Averaged Navier–Stokes (RANS) equations.

The finite volume method (FVM) is a typical numerical method to discretize
these physical equations. In FVM, the computational space is discretized into
small cells by dense grid points, separately solving every cell’s evolution at
a small timestep and finally integrating all time steps. The PDE will be
converted to a sparse linear equation at every time step, whose dimension $N$
has a linear dependency on the number of cells. In practice, sparse matrix
linear solvers such as the conjugate gradient method are available. The best
time complexity of the conjugate gradient method is $O(Ns\kappa\log
1/\epsilon)$ time complexity where $s$ is the sparsity number (the maximum
number non-zero elements in each row or column), $\kappa$ is the condition
number, and $\epsilon$ is the precision [1]. A typical problem for the FVM is
that when the problem size grows large, the computing resources will become
expensive.

Instead of using classical computers, quantum computing is a promising
computing paradigm that offers exponential acceleration over classical
computing approaches. Many quantum algorithms, including quantum
factorization[2], quantum simulation[3, 4, 5, 6], and the linear system
solvers, [7, 8, 9] have already appeared to prove this idea. Thus, we try to
accelerate a CFD solver with quantum computing. There were some works about
solving linear PDEs using a quantum computer[10, 11, 12, 13]. However, these
methods cannot be directly applied to solve CFD because the Navier–Stokes
equation is a non-linear PDE, which is not covered by these previous results.

This paper introduces a quantum solver for CFD problems (QCFD solver) based on
classical FVM. We show that with only classical input, the time cost of each
iteration step can be reduced to polylogarithmic dependency on the problem
size. This provides an exponential speedup of the FVM. The QCFD solver can
fully reproduce the result of the FVM. With the output given classically at
every time step, the QFVM is capable of the steady or unsteady problem with
similar configurations.

To apply quantum algorithms to practical problems, the conversion between the
classical and the quantum data could become a bottleneck, especially when
using the quantum linear solver (QLS) as a submodule[14]. In our algorithm,
the input and output are all classical data. To achieve this, we design a
quantum memory layout based on quantum random access memory (QRAM) [15, 16].
At the input stage, the memory layout helps to implement subprocedures
required by the QLS. At the output stage, we sample the output state and
update the memory classically. We show that these two processes, which act as
the interface between classical and quantum data, can both run in
polylogarithmic time. They enable us to integrate the quantum linear solver
submodule into the classical FVM to achieve speedup. The time complexity of
our algorithm is calculated by scaling the time between two iteration step,
which is

$O\left(\frac{(s^{3}+\log
N)s\kappa\log^{3}N}{\epsilon^{2}}\mathrm{polylog}(s\kappa/\epsilon)\right).$
(1)

Our algorithm can be compared to the classical algorithm directly. The
algorithm has classical input and output and uses the same definition of
condition number and the error threshold. For the condition number problem, we
implement a quantum version of the Jacobi preconditioner and integrate it into
the memory layout. It has the same effect as applying a Jacobi preconditioner,
which is common in the classical FVM. Therefore, both quantum and the
classical algorithm have a linear dependency (if we ignore the polylogarithmic
term) on the same condition number. We also provide evidence that the quantum
and the classical error threshold are the same under the sense of $l_{\infty}$
norm.

As a result, compared to the time complexity of the classical counterpart, the
quantum solver runs faster in terms of the problem size $N$, but slower in the
dependency on the precision with a quadratic term. Our algorithm will have
better performance when the problem size is large enough, i.e., $N\gg
1/\epsilon^{2}$. We performed a numerical simulation on the Onera M6 test case
showing that the algorithm can output correctly with such a condition.

##  II Preliminaries

###  II.1 Discretization and Linearization of the Physical Governing Equation

The typical physical governing equations (Euler, NS, RANS) have to be
linearized to apply to the FVM. In this paper, we do not focus on the detail
of the linearization. Instead, we apply the identical linearization method to
the classical algorithm and analyze the relationships among the equation
variables.

Here we take a two dimensional NS equation with compressible flow as an
example. First write down the differential form of the NS equation:

$\frac{\partial}{\partial t}\int_{\Omega}UdV+\oint_{\partial\Omega}\bm{F}\cdot
d\bm{S}=0,$ (2)

where

$\begin{array}[]{ccc}U=\left[\begin{array}[]{c}\rho\\\ \rho u\\\ \rho v\\\
\rho E\end{array}\right]&\bm{F}_{x}=\left[\begin{array}[]{c}\rho u\\\ \rho
u^{2}+p\\\ \rho uv\\\ \rho
uH\end{array}\right]&\bm{F}_{y}=\left[\begin{array}[]{c}\rho v\\\ \rho uv\\\
\rho v^{2}+p\\\ \rho vH\end{array}\right],\end{array}$ (3)

for any volumne $\Omega$ and its boundary $\partial\Omega$.

Figure 1: The grid cell around $i^{\rm th}$ point. $\bm{F}_{\partial\Omega}$
is the flux at the certain boundary; $\Delta S_{\partial\Omega}$ is the area.
$\Omega_{i}$ is the volume of this cell.

To discretize it spatially and timely, we split the space and time into small
grid cells. At the cell i and time step n, the NS equation can be discretized
to

$\frac{\Omega_{i}}{\Delta
t}\left(U^{n+1}_{i}-U^{n}_{i}\right)=-\sum_{\partial\Omega}\bm{F}_{i,\partial\Omega}^{n+1}\cdot\Delta\bm{S}_{i,\partial\Omega}^{n+1},$
(4)

where implicit Euler method is applied. We define the right hand side of the
equation (4) as the residual of this point, denoted by $R_{i}^{n+1}$. The
$\bm{F}_{i}$ is defined by the difference scheme, which is calculated by
variables $\vec{U}$ in the surrounding cells. The difference scheme gives a
relation between nodes. In this paper we define a matrix $C$ which has

$C_{i,i^{\prime}}=1$ (5)

if i and $i^{\prime}$ are related in the difference scheme. In other words,
$C_{i,i^{\prime}}=1$ means calculateing the residual at $i^{\rm th}$ node uses
the variables in $j^{\rm th}$ node. Specially, we always have $C_{i,i}=1$.

Let $\Delta\vec{U}^{n+1}=\vec{U}^{n+1}-\vec{U}^{n}$, we have

$\left(\frac{\Omega_{i}}{\Delta t}\delta_{i,i^{\prime}}+\left.\frac{\partial
R_{i,k}}{\partial U_{i^{\prime},k^{\prime}}}\right|_{U=U^{n}}\right)\Delta
U^{n}_{i^{\prime},k^{\prime}}=-R^{n}_{i,k}.$ (6)

Simply replacing $\bm{A}=\left(\frac{\Omega_{i}}{\Delta
t}\delta_{i,i^{\prime}}+\left.\frac{\partial R_{i,k}}{\partial
U_{i^{\prime},k^{\prime}}}\right|_{U=U^{n}}\right)$, we obtain a linear
equation whose solution implies the time evolution of the physical variable
$\vec{U}$.

The coefficient matrix $\bm{A}$ is a sparse matrix. From equation (4), the
$A_{i,k,i^{\prime},k^{\prime}}$ is non-zero when $i$ and $i^{\prime}$ are
related in the difference scheme ($C_{i,i^{\prime}}=1)$. The sparse number
(number of nonzero element in a row or column) is fixed by how we select the
difference scheme, denoted by $s$.

Regardless of the physical governing equation, the discretization and the
linearization following the classical FVM method do not change. We will
finally show that the spatial or time difference scheme does not affect how
this algorithm works, and only the constant-coefficient will change in the
analysis of the time complexity.

###  II.2 Quantum Algorithm with Classical Input and Output

When we use a quantum computer to cope with a practical problem, we should
always expect that the input and output are classical. Many quantum algorithms
have been proposed and claimed to be faster (exponentially or polynomially)
than their classical counterparts. However, a large portion of them only beats
classical algorithm under some theoretical limitations. A typical example is
the famous quantum linear system algorithm: Harrow-Haddism-Lloyd (HHL)
algorithm[7], which can prepare the state $|x\rangle$ encoding the solution of
the linear equation $Ax=b$. This algorithm uses $O(\log N)$ calls to linear
equation oracles, where the classical counterpart has to perform at least
$O(N)$ calls. Based on the ability to accelerate solving a linear equation,
many quantum machine learning algorithms were proposed and accelerated
exponentially over the classical algorithm. However, most of these algorithms
did not answer how to deal with real-world data or obtain a classical output.
In this paper[14], Aaronson raised a series of obstacles for applying the HHL
algorithm on quantum machine learning algorithms with real-world data. The
main problems include how to input the classical data into the quantum
computer and extract information from the output state given by the HHL
algorithm. If we hope to preserve quantum speedup, two operations are
forbidden. One is to prepare the input state $|b\rangle$ with an encoded
quantum circuit, where even reading all data entries requires $O(N)$ time; the
other is to perform sampling on the output state to extract the state to a
classical vector with $O(N)$ times measurement.

We believe the obstacles that appeared in “HHL-based” quantum machine learning
algorithms are even more challenging if we want to accelerate the FVM for CFD
problems quantumly. In the classical method, the physical governing equation
(e.g., Navier–Stokes equation) will be discretized in both time and space. At
every time step, we linearize the physical governing equation to a linear
equation, and its solution represents the evolution of the physical variables.
Finally, we can integrate all time steps to obtain a stable solution. The time
complexity of the FVM mainly depends on the time complexity of the linear
solver, which is at least $O(N)$. It is straightforward to suppose that if we
change the linear solver to a quantum version, we will have a quantum
accelerated CFD solver. However, the two problems mentioned above exist in
every step of the time integration. One is how to generate the input from the
physical variables at one time step; the other is how to update the physical
variables from the quantum linear solver’s output. In conclusion, the quantum
algorithm will fail to demonstrate quantum advantage if time complexity
requires an extra $O(N)$ multiplier.

Our proposed QFVM algorithm will consider these obstacles. We assume the input
and output of this algorithm are all classical, ensuring that this algorithm
can be run in the quantum computer without providing more input than the
classical algorithm.

###  II.3 Quantum Random Access Memory

Quantum random access memory is the storage device for the quantum computer.
As the quantum analog of RAM, QRAM allows a quantum computer to obtain
classical data with given addresses in quantum parallel. In other words, QRAM
could perform such unitary transformation:

$U_{\rm QRAM}|i\rangle_{A}|0\rangle_{D}=|i\rangle_{A}|d_{i}\rangle_{D},$ (7)

where A and D denote the address and the data registers. $d_{i}$ is a
classical data entry stored at the address $i$.

A seminal architecture called “bucket-brigade” provides an efficient way for
querying. There have been many proposed physical implementations of such
architecture, such as optical system [15], acoustics system [17], and circuit
quantum electrodynamics[18]. Our work is based on the QRAM with architecture
implemented by any of the physical systems. To eliminate the difference in
understanding the availability of the QRAM, we list all assumptions when we
apply the QRAM to our algorithm.

First, the QRAM is general to all input addresses and their superpositions,
namely $\sum c_{i}|i\rangle$. The QRAM should be an arbitrary data loader
rather than only allowing to prepare the $\sum|i\rangle|d_{i}\rangle$ state.

Second, if the address register has been prepared, performing one query costs
$O(\log N)$ time where the full data length is N.

Third, we assume a QRAM has at least a classical RAM capability, enabling
access to a single entry or overlaying it to another value with constant time.
Meanwhile, the QRAM should be compatible with a classical computer. A
classical computer can read the data in QRAM without extra cost.

Even though a real physical implementation of QRAM is hard, these assumptions
are reasonable because they do not exceed the capabilities of the previous
physical implementations.

We claim our algorithm as “classically input and output” under the sense that
the input and output of the algorithm are stored in the QRAM. Because we
believe in the compatibility of QRAM and classical computer, the problem
definition, data initialization, and post-processing of the calculation
results can all be performed in a connected classical computer.

##  III Quantum Finite Volume Method

This section introduces the quantum finite volume method (QFVM) for CFD
problems based on the implicit Euler method and classical FVM. As described in
the preliminary section (II.1), the Euler, NS, or RANS equation can be
discretized timely and spatially and finally linearized following specific
classical methods. We do not focus on these methods but only transplant them
into a quantum version based on the theory that any classical function can be
implemented in a quantum computer.

Take a two-dimensional Euler equation with compressible gas (density, X- or
Y-directional momentum and energy) as the example. The physical variable
$\vec{U}$ contains the physical properties of all grid cells. We discretize
the computing space into N grid cells, then the size of $\vec{U}$ is 4N. We
use a two-level subscript $(i,k)$ to identify a single element in $\vec{U}$ by
$U_{i,k}$, where $i$ denotes the $i^{\rm th}$ grid cell and $k$ denotes the
$k^{\rm th}$ physical characteristics.

The time evolution of $\vec{U}$ is realized by solving the linearized NS
equation using implicit Euler time iteration scheme. The calculation is
iterative, and we use a subscript n to denote the $n^{\rm th}$ iteration step.
At any time step, the equation has a general form:

$\bm{A}^{n}(\vec{U}^{n})\Delta\vec{U}^{n+1}=-\vec{R}^{n}(\vec{U}^{n}),$ (8)

where $\bm{A}$ is the Jacobian matrix and $\vec{R}$ is the residual vector.
These two parts are determined by the physical variable $\vec{U}$ at $n^{\rm
th}$ step. The unknown is $\Delta\vec{U^{n+1}}$. The system evolves one step
by calculating

$\vec{U}^{n+1}=\vec{U}^{n}+\Delta\vec{U}^{n+1}.$ (9)

Similar to the definition of $U_{i,k}$, the matrix element of $\bm{A}$ is
denoted by $A_{i,k}^{i^{\prime},k^{\prime}}$. Element in $\vec{R}$ is denoted
by $R_{i,k}$.

The initial state of the physical variable $\vec{U}^{0}$ should have been
stored in QRAM. It is given as the input of the algorithm. Another part of the
input is about the spatial grid points, which split the space into cells.
These data should also be stored in the QRAM, ready for quantum query.

Our work mainly concentrates on how to bridge the gap between quantum
processes and classical processes. On top of our work, we design a QRAM-based
memory layout inspired by previous works[19]. The details about the memory
layout are introduced in section IV.

Alike the classical program, the QRAM stores the physical variable
$\vec{U}^{n}$ to construct the linear equation at the step n. Besides this, we
also prepare the residual vector $\vec{R}$ and a sum tree With the memory
layout design, one can construct three quantum subprocedures $O_{A}$, $O_{b}$,
and $O_{l}$ required by the QLS. They encode the linear equation as unitary
transforms, that is:

$O_{A}|i,k,i^{\prime},k^{\prime}\rangle=|i,k,i^{\prime},k^{\prime}\rangle|A_{i,k}^{i^{\prime},k^{\prime}}\rangle,$
(10)

which encodes the Jacobian matrix’s element, and

$O_{b}|i,k\rangle=|i,k\rangle|R_{i,k}\rangle,$ (11)

which encodes the residual vector’s element, and

$O_{l}|i,p\rangle=|i,C_{i}(p)\rangle,$ (12)

which encodes the $p^{\rm th}$ related cell in the difference scheme.

$O_{l}$ can be implemented directly by the geometry definition input with
constant queries. $O_{A}$ also requires querying the geometry definition to
obtain the $O(s)$ number of $U_{i}$ at the related cells. Implementation of
$O_{b}$ is introduced in section IV.2.

With these three quantum subprocedures, QLS outputs a solution
$|\vec{u}\rangle=|\frac{\Delta\vec{U}^{n+1}}{\lVert\Delta\vec{U}^{n+1}\rVert_{2}}\rangle$,
a normalized solution of the linear equation. Now we are able to construct a
procedure $P$ which can prepare the $|\vec{u}$ with sublinear time.

Taking $P$ as the input of the $l_{\infty}$ tomography algorithm, we can
obtain a classical vector $\tilde{u}$ which is $\epsilon$-close to the quantum
solution $\vec{u}$. This algorithm requires to run $P$ and its controlled
version by $O(\frac{\log N}{\epsilon^{2}})$ many times. As a result, obtaining
$\tilde{u}$ requires sublinear time.

Meanwhile, we can apply amplitude estimation[20] to $P$ to obtain the
normalization factor $c_{l}=\lVert\Delta\vec{U}^{n+1}\rVert_{2}$. we use
$\tilde{u}$ and $c_{l}$ to update the sum tree and finish this iteration step.

For a steady problem, the computing stops in two cases. One is when the
residual is smaller than the convergence limit $\epsilon$, which can be
extracted from the top of the tree. Another is when reaching the maximum
iteration steps. After stopping, the output of this algorithm is stored in the
physical variable area.

##  IV Quantum memory Layout for QFVM Algorithm

To reduce data transfer costs between quantum and classical data, we design a
memory layout to efficiently prepare the residual vector state and prepare the
oracular input of the quantum linear solver.

Figure 2: Schematic of the quantum memory layout. (a)Three parts of the QRAM
memory layout: geometry definition area, which holds the input of the problem;
physical variable area holding the $\vec{U}^{n}$; and the residual sum tree.
(b)The linear structure of the geometry definition area. This area is formed
with N blocks. The $i^{\rm th}$ block holds s related indices where each of
$i^{\prime}=i_{k}$ ($0\leqslant k<s$) satisfies $C_{i,i^{\prime}}=1$. (c)The
linear structure of the physical variable area. Each block the physical
variables $\vec{U}_{i}$ at the cell i. (d)The binary tree structure of the
residual sum tree. The tree’s leaves are the components of the residual vector
$\vec{R}^{n}$. Then for each level, we sum up the square of every two nodes.
The tree root is $\lVert\vec{R}\rVert^{2}$.

The QRAM stores three kinds of data. The first is the geometry data of the
problem definition. From this, we can query the connection relation between
cells (geometry definition) in quantum parallel. These are constant during the
calculation and used for constructing the linear equation. Querying this part
is equivalent to this unitary transform:

$\mathcal{P}_{s}|i\rangle|j\rangle=|i\rangle|i_{j}\rangle.$ (13)

The second is the vector of the physical variable at the step n, namely
$\vec{U}^{n}$. From this part we can implement

$\mathcal{P}_{U}|i,k\rangle|0\rangle=|i,k\rangle|U_{i,k}\rangle.$ (14)

The third is the residual vector $\vec{R}^{n}$ and its pre-computed sum. We
call it the sum tree. From the bottom to the top of the binary tree, every
node stored an integration necessary to prepare the residual state. We define
each node’s address as $a_{r}(p)$ from the top of the sum tree. $p$ is a
binary string where every digit represents the left/right branch with 0/1. For
example, $a_{r}(0)$ is the address of the left child of the root; $a_{r}(0,1)$
is the right child of the node at $a_{r}(0)$. Specially, we directly use
$a_{r}$ to represent the root’s address. We place all nodes in the QRAM
linearly, where each node’s address $a_{r}(p)$ can be computed efficiently.
The data contained in address $a_{r}(p)$ is denoted by $S_{R}(p)$. Thus we can
perform such unitary transform:

$\mathcal{P}_{R}|a_{r}(p)\rangle|0\rangle=|a_{r}(p)\rangle|S_{R}(p)\rangle.$
(15)

The second and the third part will be initialized once at the beginning of the
calculation. They will be continuously updated through the calculation. The
final steady result, which is $\vec{U}$ at the last iteration step, is also
stored in them.

The diagram of the memory layout is demonstrated in figure 2. The functions of
this memory layout are as follows:

  * •

Initialize the sum tree will classically access the QRAM $O(N)$ times;

  * •

With access to the sum tree, one can prepare $|R\rangle$ with $O(\log N)$
times of quantum queries to the QRAM;

  * •

If the sparsity of the Jacobian matrix is $s$, updating a single entry of
$\vec{U}$ will classical access the QRAM $O(s\log N)$ times to update the sum
tree;

  * •

The normalization factor between the residual vector and its corresponding
quantum state can be obtained with one quantum query;

  * •

Evolution from $U^{n}$ to $U^{n+1}$ will cost $O(C\log^{2}N/\epsilon^{2})$
time.

We clearly state the difference between “classical access” and “quantum
query”. Classical accessing means only one data entry is read and modify at a
time. Quantum query means a unitary transform is performed, simultaneously
extracting many data queries into one quantum register using the superposition
addresses.

In the following part, we will show how these features are realized.

###  IV.1 Initialization

The initialization process fills the QRAM as the memory layout scheme shows.
The first step is to write in the initial physical variable $\vec{U}$ and fill
the tree with the wanted sum.

The initialization is entirely a classical process. Along with writing the
data, we should also record the memory layout to quickly obtain the memory
address of every data entry in constant time.

There may be some concentrations about whether such $O(N)$ preparation time
will cause the vanishment of this algorithm’s speedup. We believe after
considering the time consumption of initialization, the quantum speedup still
preserves. We analyze the three things that contribute to initialization time.

One is the calculation of the residual vector $\vec{R}$. The calculation of
this vector exists at every step of the classical FVM in a CFD problem. Even
in classical algorithms, this part is not the bottleneck of the time. Our
algorithm only calculates the residual once initially, which will consume much
less than a classical algorithm does.

The second is the fill of the sum tree. To fill a sum tree only requires
repeatedly adding the sum of the square of the residual vector. It is natural
to think this process is easier than the calculation of the residual vector.

The third is about the cost of accessing the QRAM classically. As we have
mentioned in section II.3, we assume the QRAM has the near capability of RAM,
which allows the access to be performed in constant time.

###  IV.2 Preparation of the residual vector state $|\vec{R}\rangle$

According to the method described in [21], the state $|\vec{R}\rangle$ could
be prepared efficiently because we have access to all wanted sums of the
vector. We pre-compute them in the sum tree, so the preparation can be
realized by querying the sum tree $\mathcal{P}_{R}$.

The first step is to query the tree root and its left child node, then
calculate the rotating angle at this step:

$|a_{r}(0)\rangle|S_{R}(0)\rangle|a_{r}\rangle|S_{R}\rangle|\theta\rangle|0\rangle,$
(16)

where $S_{R}(0)=\sum_{i\in[0..N/2-1]}R_{i}^{2}$,
$S_{R}=\sum_{i\in[0..N-1]}R_{i}^{2}$, $\theta=\arccos\frac{S_{R}(0)}{S_{R}}$.

Now perform a conditional rotation and uncompute, we have:

$\cos\theta|0\rangle+\sin\theta|1\rangle.$ (17)

Add another qubit, perform Hadamard gate on it, we have

$\displaystyle(\cos\theta|0\rangle+\sin\theta|1\rangle)\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$
(18)
$\displaystyle=\frac{1}{\sqrt{2}}(\cos\theta|00\rangle+\sin\theta|10\rangle)$
(19)

Then we iteratively perform the query, computing the rotating angles and
conditional rotation. At $k^{\rm th}$ step, we have the state

$\sum c_{i}^{k}|i\rangle\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle).$ (20)

Computing the addresses of $|i,0\rangle$, $|i,1\rangle$ , we obtain the real
addresses $a_{r}(i,0)$ and $a_{r}(i,1)$ in the QRAM, that is

$\sum
c_{i}^{k}|i\rangle\frac{1}{\sqrt{2}}(|0\rangle|a_{r}(i,0)\rangle+|1\rangle|a_{r}(i,1)\rangle,$
(21)

then query to the $\mathcal{P}_{R}$ to obtain the rotating angles and
uncompute extra registers. Finally, after performing conditional rotation, we
step to

$\sum c_{i}^{k+1}|i\rangle$ (22)

Repeatedly performing this process, we can efficiently prepare the residual
state $|\vec{R}\rangle=-\sum
R_{i}|i\rangle/\sqrt{\sum_{j\in[0..N-1]}R_{j}^{2}}$ with the help of the sum
tree.

###  IV.3 Updating a single entry of the physical variable

From step n to n+1, the physical variable $\vec{U}$ should be updated. Same as
the computing $\vec{R}$ and its sum from $\vec{U}$, only the residual on the
related cells would change. From the tree leaves, we change the all residual
$R_{i^{\prime},k^{\prime}}$ related to the changed $U_{i,k}$ with
$C_{i,i^{\prime}}=1$. After these residual vector entries change, we again
compute the sum tree from the leaves to the root and update correspondingly.
The number of updated nodes will not exceed the number of the multiplication
of the related residual entries $O(s)$ and the number of layers of the sum
tree $\log N$.

As a result, the cost of update one entry of $\vec{U}$ is less than $O(s\log
N)$.

###  IV.4 Sampling the solution state and update the QRAM

The QLS outputs the solution as a quantum state. However, we cannot directly
extract it to a classical vector to update $\vec{U}$. We need to cope with two
problems: first is to decide the normalization factor of the solution; the
second is to convert the quantum state to a classical vector.

The tomography algorithm only produces a normalized vector
$\Delta\tilde{U}^{n+1}$. We should also obtain all the normalized factor in
the algorithm to get a real update vector of $\vec{U}^{n+1}$. QLS produces two
factors. First is $c_{b}$, which is generated when preparing $|b\rangle$ and
can be obtained from the data structure described above. The second is $c_{l}$
which is generated from the QLS, because matrix inversion is usually not
unitary and the raw solution $A^{-1}b$ is not normalized. With amplitude
estimation, we can compute the probability $p_{l}$ and then obtain the factor
by $c_{l}=\alpha\sqrt{p_{l}}$, where $\alpha$ is a constant in the QLS.
Obtaining the normalization factors will not affect the asymptotic time
complexity of the algorithm.

Combining these two factors $c=c_{b}c_{l}$, we obtain the norm of solution
$\lVert\Delta\vec{U}\rVert_{2}$, which implies the variation updated on the
target vector $\vec{U}$ in the CFD solver. When the norm of the variation is
smaller than a preset threshold, we can stop the iteration and return the
result.

With $l_{\infty}$ tomography[22] we could efficiently produce an
$l_{\infty}$-close sample $\Delta\tilde{U}^{n+1}$ of a real-valued quantum
state $|\Delta\vec{U}^{n+1}\rangle$ with $O(\frac{\log N}{\epsilon^{2}})$.
This sampling algorithm has a logarithmic dependency on N, enabling each
iteration step of our algorithm to run at polylogarithmic time complexity over
the input size N.

Updating the QRAM from the sampled vector is also efficient. The $l_{\infty}$
tomography algorithm produces a sparse classical vector with not more than
$O(\log N/\epsilon^{2})$ non-zero elements, which means the updating will be
performed for less than $O(\log N/\epsilon^{2})$ times to update the QRAM
$\mathcal{P}_{U}$ in one iteration step. This results in efficiency in both
sampling and updating.

##  V Implementation of Quantum Preconditioner

The condition number of the linear equation represents to what extent the
solution can be affected by the pertubation on the right-hand-side vector. The
condition number is defined as:

$\kappa(\bm{A})=\frac{|\lambda_{max}|}{|\lambda_{min}|},$ (23)

where $|\lambda_{max}|$ and $|\lambda_{min}|$ is the maximum/minimum absolute
of the eigenvalues of $\bm{A}$. When the condition number is large, we say the
equation is ill-conditioned, requiring high precision and time complexity to
solve. The time complexity of the classical sparse linear solver has a
dependency on the condition number. For example, the time complexity of the
conjugate gradient method is $O(\kappa sN\log 1/\epsilon)$. The QLS used in
our algorithm also has a linear dependency on the condition number.

Preconditioner is a pre-processing method that can reduce the condition number
of the equation. If we have a matrix P such that
$\kappa(\bm{PA})<\kappa(\bm{A})$, we can transform this equation as:

$\bm{A}\vec{x}=\vec{b}\Leftrightarrow\bm{PA}\vec{x}=\bm{P}\vec{b}.$ (24)

Preconditioners are constructed from the raw equation, and there have been
many types of preconditioners. However, not all classical preconditioners
could be directly transplanted to quantum versions. First, the matrix
multiplication by the preconditioner should be computed efficiently, namely
within $O(\mathrm{polylog}(N))$ time. Second, the preconditioned matrix should
also be sparse; otherwise, it cannot be efficiently solved by the QLS. Some
preconditioners suitable for QLS have already proposed in [23] and [24]. We
here display an example preconditioner: the blockwise Jacobi preconditioner,
which is widely used in the classical CFD solver. We implement the blockwise
Jacobi preconditioner in our algorithm without affecting the asymptotic
complexity on the problem size N.

#### Apply Jacobi preconditioner to subprocedures

Jacobi preconditioner uses the inverse of the diagonal block. For the raw
linear equation, we construct subprocedures as the input of the QLS. The
preconditioned equation has a different matrix and vector; therefore, these
subprocedures should be modified.

Let $\tilde{\bm{A}}=\bm{PA}$ and $\vec{R^{\prime}}=\bm{P}\vec{R}$, where
$\bm{P}$ is the Jacobi preconditioner of the matrix $\bm{A}$. The element of
the $\bm{P}$ is

$\bm{P}_{i,k}^{i^{\prime},k^{\prime}}=\delta_{i,i^{\prime}}B_{i,k}^{i^{\prime},k^{\prime}},$
(25)

where $\bm{B}_{i,\\_}^{i^{\prime},\\_}$ represents the inverse of the block
$\bm{A}_{i,\\_}^{i^{\prime},\\_}$.

The element of $\tilde{\bm{A}}$ is

$\tilde{A}_{i,k}^{i^{\prime},k^{\prime}}=\sum_{j\in[0..N-1];l\in[0..n_{\rm
var}-1]}\bm{P}_{j,l}^{i^{\prime},k^{\prime}}\bm{A}_{i^{\prime},k^{\prime}}^{j,l}.$
(26)

The Jacobi preconditioner $\bm{P}$ is blockwise diagonal. We can simplify the
equation (26) as:

$\tilde{A}_{i,k}^{i^{\prime},k^{\prime}}=\sum_{l\in[0..n_{\rm
var}-1}]\bm{P}_{i,l}^{i,k^{\prime}}\bm{A}_{i^{\prime},k^{\prime}}^{i,l}.$ (27)

This implies that computing a single element of $\tilde{\bm{A}}$ requires to
queries $n_{\rm var}$ elements of $\bm{A}$. Another fact is that the sparsity
matrices of $\tilde{\bm{A}}$ and $\bm{A}$ are the same when they are symmetric
to the diagonal line. This is often true because in the difference scheme, $i$
and $i^{\prime}$ are related so that $A_{i,k}^{i^{\prime},k^{\prime}}$ and
$A_{i^{\prime},k}^{i,k^{\prime}}$ are all non-zero elements.

When it is efficient to implement $O_{A}$, $O_{A}^{\prime}$ will also be
efficient to implement. That is

$O_{A}^{\prime}|i,k,i^{\prime},k^{\prime}\rangle|0\rangle=|i,k,i^{\prime},k^{\prime}\rangle|\tilde{A}_{i,k}^{i^{\prime},k^{\prime}}\rangle.$
(28)

First we query all elements required for computing the inverse at
$(i,i^{\prime})$ block and compute the inverse, we have

$|i\rangle|B_{i,\\_}^{i,\\_}\rangle_{\rm regs(A)}.$ (29)

We use $B_{i,\\_}^{i^{\prime},\\_}$ to represent a matrix block with $n_{\rm
var}^{2}$ elements. The subscript “regs(A)” mean we require a group of quantum
registers to hold this matrix, marked by A.

The corresponding block in $\bm{A}$ is also queried,

$|i,i^{\prime}\rangle|A_{i,\\_}^{i^{\prime},\\_}\rangle_{\rm regs(B)},$ (30)

Combining two register group A and B, we obtain the wanted element
$A_{i,k,i^{\prime},k^{\prime}}$. Computing one element requires $n_{\rm
var}^{2}$ times of calls to the $O_{A}$, namely propotional to $O(s^{2})$.

From the above derivation, the sparsity of $\tilde{\bm{A}}$ is same with the
$\bm{A}$. Therefore $O_{l}$ remains unchanged.

The subprocedure $O_{b}$ should also be modified. In the original description
of the memory layout, the sum tree stores the pre-computed residual vector. In
the preconditioned version, the $\vec{R}$ is replaced by $\vec{R^{\prime}}$.
When the $\vec{U}$ changes according to the sampling results, we need to
compute the preconditioned residual vector and update the sum tree. Computing
any element of $R^{\prime}_{i,k}$ is still related to all connected cells,
which is

$R^{\prime}_{i,k}=\sum_{k^{\prime}\in[0..n_{\rm
var}-1]}\bm{P}_{i,k^{\prime}}^{i,k}R_{i,k^{\prime}}$ (31)

The complexity of this process is also contributed by computing the diagonal
block’s inverse of $\bm{A}$.

From the modified sum tree, constructing such $O_{b}^{\prime}$ is similar to
$O_{A}^{\prime}$.

##  VI Run Time Analysis

The time cost for the QCFD algorithm has two main contributions. One is the
cost of the initialization of the memory layout (initialization cost); the
other is the time complexity between two iteration steps (evolution cost).

#### Initialization cost

In section IV.1, we estimate the time cost of the initialization stage. The
conclusion is that the initialization stage has $O(N)$ time complexity.
However, under the assumption about the QRAM’s capability (see section II.3),
the initialization cost would not cost much more than the preprocessing stage
of the classical FVM. Because the initialization of QCFD only processes once,
but classical FVM has to preprocess it as many times as the iteration steps,
we believe this cost would not become the bottleneck of the QFVM algorithm.

#### Evolution cost

Every evolution stage. In [9], Childs et al. provide a linear solver algorithm
with logarithmic dependence on precision. They show that the query complexity
of this algorithm of $O_{A}$, $O_{l}$ and $O_{b}$ are
$O\left(s\kappa\mathrm{polylog}(\frac{s\kappa}{\epsilon})\right)$. Now we
start to analyze the time complexity of constructing these subprograms from
the initial problem settings.

According to the results in the previous sections, the number of queries to
QRAM for implementing $O_{A}$, $O_{l}$, and $O_{b}$ is $O(s)$, $O(1)$, and
$O(\log N)$, correspondingly. The time complexity of preconditioned
$O_{A}^{\prime}$ has a multiplier of $O(s^{3})$ contributed by computing the
inverse of the diagonal blocks of $\bm{(}A)$. The preconditioned
$O_{b}^{\prime}$ has the same complexity as $O_{b}$.

Now consider the time cost of sampling and updating. We run the QLS with
$O(\frac{\log N}{\epsilon^{2}})$ times to obtain an $l_{\infty}$-close
classical vector. This becomes another multiplier to the time complexity of
the quantum procedure.

The last multiplier is the cost of querying the QRAM. As we have assumed, the
QRAM use $O(\log N)$ time to perform one query. By composing these results,
the time complexity of the quantum procedure is

$O\left(\frac{(s^{3}+\log
N)s\kappa\log^{3}N}{\epsilon^{2}}\mathrm{polylog}(s\kappa/\epsilon)\right).$
(32)

The final step is to update the sum tree. Updating a preconditioned residual
tree has two steps. One is to compute the preconditioned residual, where each
term will involve in another inversion of the matrix $A$, which is $O(s^{3})$;
the other is to update the tree from bottom to the top, this involves $O(\log
N)$ times for one change in the bottom of the tree. While at most
$O(\frac{\log N}{\epsilon^{2}})$ terms of $\vec{U}$ changes, the time cost of
updating the tree is $O(s^{3}\log^{2}N/\epsilon^{2})$. The total time
complexity is the addition of the quantum and the classical procedure. Because
the quantum procedure’s complexity is asymptotically greater than the
classical’s, we conclude that the evolution time cost has the time complexity
shown in equation (32).

The classical counterpart’s time complexity is $O(Ns\kappa\log 1/\epsilon)$
when using CG as the linear solver. Our algorithm outperforms the classical
algorithm on the problem size’s dependency but has worse performance when the
problem requires high precision. When the problem size $N$ and the requirement
of the precision $\epsilon$ has such relation $N\gg 1/\epsilon^{2}$, the
quantum algorithm will potentially have better performance on time.

##  VII Numerical Simulation

To prove our algorithm’s effectiveness, we performed simulation on a test case
to find whether a big size problem with low error sensitivity exists.

An open-source classical CFD simulation software, SU2[25], is selected as the
example. We appended quantum error in the simulation process, output the
evolution history, and compared it to the classical solver. The error is
implemented by biasing the solution vector with a randomized error to emulate
the sampling process of the $l_{\infty}$ tomography. Except for the error, we
preserved all physical problem configurations, including the mesh setting and
physical parameters like temperature or Mach number. The multigrid option was
turned off so that the linear equation will only be solved once in every
iteration step. The linear solver parameters are chosen to be as precise as
possible to emulate the case that the quantum linear solver produces the
result accurately.

We select the inviscid flow around Onera M6 airfoil as the test case, which
has 108396 grid points, and the physical governing equation is the three-
dimensional Euler equation. We compare the classical result with different
error settings: from 5e-2 to 1e-5. The results are displayed in figure 2. When
the error is set to 5e-2 quickly diverges. Except that, all cases converge
correctly. The black line is the classical baseline. The maximum convergence
error $\epsilon_{1}$ of this case is between 5e-2 and 3e-2, where quantum
advantage preserves.

To simulate the quantum effect, we biased the solution vector with a specific
error to emulate the sampling process of the $l_{\infty}$ tomography. This is
implemented as a postprocessor of the linear solver.

Figure 3: Results of the simulation. The test case is the three-dimensional
inviscid flow around the Onera M6 airfoil. a) Convergence history of the test
case of the quantum solver with different error setting, compared to the
classical solver as the baseline (use BCGStab as the linear solver). The error
is set from 5e-2 to 1e-5. Except $\epsilon$=5e-2, all other cases converge
correctly. The maximum convergence error is approximately 3e-2. b) Flow
density around the airfoil at quantum error $\epsilon$=3e-2. This result is
the same as the classical solver’s result displayed in subfigure d). c) Bad
solution at $\epsilon$=3e-2. The flow density at 20 iteration steps is not
correctly computed. d) The solution output by the classical solver.

We focus on the convergence history of test cases and define two kinds of
error thresholds. The first is “maximum convergence error $\epsilon_{1}$”, the
maximal value which allows the convergence. The problem will quickly diverge
if the error is larger than the first threshold. The second is “maximum stable
error $\epsilon_{2}$”, where the problem will have the same convergence
history if the error is less than this threshold. The convergence speed will
be gradually slow when the error gets larger between the first and the second
threshold. If the maximum stable error is large, we define this kind of
problem as “quantum friendly” because the quantum algorithm is much likely to
run faster than the classical. In this case, the problem size is large,
satisfying $N\gg 1/\epsilon^{2}$. This provides evidence that quantum
advantage can be realized in CFD problems.

##  VIII Error Analysis

The time complexity of the QFVM has better performance on the number of grid
cells $N$ but worse on the precision $\epsilon$, which implies that the
problem size should be large enough to show the quantum advantage. On the
other side, the numerical simulation shows that the precision should be small.
Otherwise, the time integration will not converge. There is the problem: if
the precision requirement has some dependency on the problem size, the quantum
acceleration will decrease or even vanish. In this section, we will provide
evidence that the precision will not grow with the problem size.

First, we calculate the total error generated by the quantum sampling with
error bound $\epsilon$ specified. At one step, we define the physical variable
$\vec{U}$ and its update $\Delta\vec{U}$. In QFVM, the quantum process outputs
a quantum state $|u\rangle$ which is propotional to $\Delta\vec{U}$,

$\Delta\vec{U}=\lVert\Delta\vec{U}\rVert_{2}u.$ (33)

The $l_{\infty}$ tomography outputs a classical vector $\tilde{u}$ which is
$\epsilon$-close to $u$. At any index $i$, we have

$u_{i}=\tilde{u}_{i}+e_{i},$ (34)

where the error term $|e_{i}|<\epsilon$.

Now we consider the amplitude of the $e_{i}$. When performing $l_{\infty}$
tomography, the output vector is a sample from the multinomial distribution
where the sampling number $M=C\log N/\epsilon^{2}$ and the probability
distribution ($|u_{0}|^{2}$, $|u_{1}|^{2}$, … $|u_{N-1}|^{2}$). At any term,
the standard error of such sample is:
$\sigma_{i}=\sqrt{M|u_{i}|^{2}(1-|u_{i}|^{2})}$. When $|u_{i}|$ is small
enough, we have $\sigma_{i}\sim|u_{i}|\sqrt{N}$. Now we assume the error
$e_{i}$ is approximatedly linear dependent on the standard error $\sigma_{i}$,
thus we have

$e_{i}=O(\sigma_{i})\sim O(|u_{i}|\sqrt{N}).$ (35)

The update vector output by the QFVM should be multiplied by
$\lVert\Delta\vec{U}\rVert_{2}$. As a result, the total error will be
amplified by this coefficient.

$E_{i}=\lVert\Delta U\rVert_{2}e_{i}.$ (36)

Compare two cases describing the same problem where one has $N$ cells and the
other has $kN$ (mark the variables with extra prime, e.g. $u^{\prime}$). We
can assume the distribution of $\Delta\vec{U}$ and $\Delta\vec{U}^{\prime}$ is
the same because the physical characteristic does not change. From this, we
have

$\lVert\Delta U^{\prime}\rVert_{2}^{2}=k\lVert\Delta U\rVert_{2}^{2},$ (37)

because only the vector size changes to $k$ times. From the definition of $u$
(equation (33)), this results in the decrease of the amplitude of the $u$,
i.e.

$u_{i}^{\prime}=\frac{1}{\sqrt{k}}u_{i}.$ (38)

Combining equation (35), (36) and (38), we obtain that $E_{i}=O(U_{i})$. This
result implies that the total error generated by the quantum sampling will not
change over the problem size $N$.

##  IX Conclusion

This paper developed a quantum version of the finite volume method (QFVM) for
solving the CFD problem with substantial speedup. Identical to the classical
FVM, the QFVM is iterative and can output the evolution history of the fluid
flow in the computing space. The input and output of the algorithm are both
classical data. In conclusion, we give the time complexity of between each
step is $O\left(\frac{(s^{3}+\log
N)s\kappa\log^{3}N}{\epsilon^{2}}\mathrm{polylog}(s\kappa/\epsilon)\right)$,
where $N$ is the number of cells (the problem size), $\kappa$ is the condition
number of the linear equation, $s$ is the sparsity of the Jacobian matrix, and
$\epsilon$ is the precision threshold of the output. Compared to the classical
solver’s best time complexity, which has a linear dependency on the grid size
$N$, this algorithm is exponentially faster. The speedup would be significant
when we choose an extremely large $N$, which allows the quantum computer to
solve complex CFD problems, such as tackling a larger space or finer mesh.

We use the QLS to accelerate the solution of the linear equation. However,
previous works often ignore the cost of transferring data between quantum and
classical computers. To achieve an efficient transfer, we design a memory
layout based on the QRAM. The memory layout stores the problem configuration
and the internal result. With a quantum parallel query, one can implement
quantum subprocedures necessary for the QLS; it is also updated efficiently
when the algorithm iterates over steps.

Numerical simulations are conducted to check whether the final result is
affected by the quantum error. The problem size of the test case is $1e5$, and
the QCFD solver converges correctly at $\epsilon$=3e-2, after around 200
steps. This result shows that CFD could have strong error tolerance with large
problem size; therefore, the quantum advantage preserves.

Our future work will focus on how precision affects the final result and the
convergence history, and how to optimize it. We believe that the quantum
computer will show its advantage in solving a more complex CFD problem
shortly.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(Grants Nos. 11625419), the National Key Research and Development Program of
China (Grant No. 2016YFA0301700), the Strategic Priority Research Program of
the Chinese Academy of Sciences (Grant No. XDB24030600), and the Anhui
Initiative in Quantum Information Technologies (Grants No. AHY0800000).

## References

  * Shewchuk [1994] J. R. Shewchuk,  _An Introduction to the Conjugate Gradient Method Without the Agonizing Pain_ (Carnegie Mellon University, 1994).
  * Shor [1999] P. W. Shor, Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer, SIAM review 41, 303 (1999).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Reviews of Modern Physics 86, 153 (2014).
  * Berry  _et al._ [2014] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Exponential improvement in precision for simulating sparse hamiltonians, in  _Proceedings of the Forty-sixth Annual ACM Symposium on Theory of Computing_, STOC ’14 (ACM, New York, NY, USA, 2014) pp. 283–292.
  * Berry  _et al._ [2015] D. W. Berry, A. M. Childs, and R. Kothari, Hamiltonian simulation with nearly optimal dependence on all parameters, in  _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ (IEEE, 2015) pp. 792–809.
  * O’Malley  _et al._ [2016] P. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, _et al._ , Scalable quantum simulation of molecular energies, Physical Review X 6, 031007 (2016).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Physical review letters 103, 150502 (2009).
  * Ambainis [2010] A. Ambainis, Variable time amplitude amplification and a faster quantum algorithm for solving systems of linear equations, arXiv preprint arXiv:1010.4458 (2010).
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46, 1920 (2017), https://doi.org/10.1137/16M1087072 .
  * Childs  _et al._ [2020] A. M. Childs, J.-P. Liu, and A. Ostrander, High-precision quantum algorithms for partial differential equations, arXiv preprint arXiv:2002.07868 (2020).
  * Fillion-Gourdeau and Lorin [2017] F. Fillion-Gourdeau and E. Lorin, Simple digital quantum algorithm for symmetric first order linear hyperbolic systems, arXiv preprint arXiv:1705.09361 (2017).
  * Costa  _et al._ [2017] P. C. S. Costa, S. Jordan, and A. Ostrander, Quantum algorithm for simulating the wave equation, arXiv preprint arXiv:1711.05394 10.1103/PhysRevA.99.012323 (2017).
  * Cao  _et al._ [2012] Y. Cao, A. Papageorgiou, I. Petras, J. Traub, and S. Kais, Quantum algorithm and circuit design solving the poisson equation, arXiv preprint arXiv:1207.2485 10.1088/1367-2630/15/1/013021 (2012).
  * Aaronson [2015] S. Aaronson, Read the fine print, Nature Physics 11, 291 (2015).
  * Giovannetti  _et al._ [2008a] V. Giovannetti, S. Lloyd, and L. Maccone, Architectures for a quantum random access memory, Physical Review A 78, 052310 (2008a).
  * Giovannetti  _et al._ [2008b] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum random access memory, Physical review letters 100, 160501 (2008b).
  * Hann  _et al._ [2019] C. T. Hann, C.-L. Zou, Y. Zhang, Y. Chu, R. J. Schoelkopf, S. M. Girvin, and L. Jiang, Hardware-efficient quantum random access memory with hybrid quantum acoustic systems, arXiv preprint arXiv:1906.11340 (2019).
  * Naik  _et al._ [2017] R. Naik, N. Leung, S. Chakram, P. Groszkowski, Y. Lu, N. Earnest, D. McKay, J. Koch, and D. Schuster, Random access quantum information processors using multimode circuit quantum electrodynamics, Nature communications 8, 1904 (2017).
  * Kerenedis and Prakash [2016] I. Kerenedis and A. Prakash, Quantum recommendation systems, (2016).
  * Brassard  _et al._ [2000] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, Quantum amplitude amplification and estimation, (2000), arXiv:quant-ph/0005055 .
  * Grover and Rudolph [2002] L. Grover and T. Rudolph, Creating superpositions that correspond to efficiently integrable probability distributions, (2002).
  * Iordanis  _et al._ [2019] K. Iordanis, L. Jonas, and A. Prakash, quantum algorithms for deep convolutional neural networks, arXiv preprint arXiv:1911.01117 (2019).
  * Clader  _et al._ [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Preconditioned quantum linear system algorithm 10.1103/PhysRevLett.110.250504 (2013), arXiv:1301.2340 .
  * Shao and Xiang [2018] C. Shao and H. Xiang, Quantum circulant preconditioner for linear system of equations 10.1103/PhysRevA.98.062321 (2018), arXiv:1807.04563 .
  * Economon  _et al._ [2016] T. D. Economon, F. Palacios, S. R. Copeland, T. W. Lukaczyk, and J. J. Alonso, SU2: An Open-Source Suite for Multiphysics Simulation and Design, AIAA Journal 54, 828 (2016).