2.1 One-dimensional Vlasov–Ampère system

Electrostatic oscillations of a one-dimensional electron plasma can be described by the Vlasov–Ampère system

(2.1a)\begin{gather} \partial_t f + v\partial_x f - \frac{e}{m_e}E\,\partial_{v}f = \mathcal{F}, \end{gather}

(2.1b)\begin{gather}\partial_t E - 4{\rm \pi} e \int v f \,{\rm d} v ={-} 4{\rm \pi} j^{(S)}, \end{gather}

where $E(t, x)$ is the electric field, $f(t, x, v)$ is the electron probability distribution, $t$ is time, $x$ is the coordinate in the physical space, $v$ is the coordinate in the velocity space, $e > 0$ is the elementary charge and $m_e$ is the electron mass. We have also introduced a fixed source term $\mathcal {F} = \mathcal {F}(x, v)$ (balanced by particle losses through the plasma boundaries), which is explained below. The term $j^{(S)}$ represents a prescribed source current that drives plasma oscillations. The corresponding source charge density $\rho ^{(S)}$ can be inferred from $j^{(S)}$ using the charge conservation law

(2.2)\begin{equation} \partial_t \rho^{(S)} + \partial_xj^{(S)} = 0. \end{equation}

Let us split the electron distribution into the background distribution $F$ and a perturbation $g$:

(2.3)\begin{equation} f(t,x,v) = F(x,v) + g(t,x,v). \end{equation}

We assume that $g$ is small, so the system can be linearized in $g$. Also, because we assume a neutral plasma, the system does not have a background electric field, so the stationary background distribution satisfies $v\partial _x F = \mathcal {F}$. Provided that $\mathcal {F}$ depends on $x$, $F$ can be spatially inhomogeneous. At the same time, since $\mathcal {F}$ is fixed, it does not enter the equation for $g$. This leads to the following linearized equations:Footnote ²

(2.4a)\begin{gather} \partial_t g + v\partial_xg - E\partial_v F = 0, \end{gather}

(2.4b)\begin{gather}\partial_t E - \int v g \,{\rm d} v ={-} j^{(S)}. \end{gather}

Here, $v$ is normalized to $v_{\rm th}$, $x$ is normalized to the electron Debye length $\lambda _D$ and the time is normalized to the inverse electron plasma frequency $\omega _{p}^{-1}$, where

(2.5a–c)\begin{equation} \omega_{p} = \sqrt{\frac{4{\rm \pi} e^2n_{\rm ref} }{m_e}},\quad v_{\rm th} = \sqrt{\frac{T_{\rm ref}}{m_e}},\quad \lambda_{D} = \frac{v_{\rm th}}{\omega_{p}}, \end{equation}

where $n_{\rm ref}$ and $T_{\rm ref}$ are some fixed values of the electron density and temperature, respectively. The distribution functions $g$ and $F$ are normalized to $n_{\rm ref}/v_{\rm th}$, and $E$ is normalized to $T_{\rm ref}/(e\lambda _D)$. We also assume that $F$ is Maxwellian, i.e.

(2.6)\begin{equation} F(x,v) = \frac{n(x)}{\sqrt{2{\rm \pi} T(x)}} \exp\left(-\frac{v^2}{2T(x)}\right), \end{equation}

where the background density $n$ and temperature $T$ are normalized to $n_{\rm ref}$ and $T_{\rm ref}$, respectively. An analytical description of the system (2.4) is presented in Appendix A.

2.2 Boundary-value problem

To reformulate (2.4) as a boundary-value problem, we consider a source oscillating at a constant real frequency $\omega _0$:

(2.7)\begin{equation} j^{(S)}(t,x) = j^{(S)}(x) \exp({-\mathrm{i}\omega_0 t}). \end{equation}

Assuming $g, E \propto \exp (-\mathrm {i}\omega _0 t)$, one can recast (2.4) as

(2.8a)\begin{gather} \mathrm{i}\omega_0 g - v\partial_x g - v H E = 0, \end{gather}

(2.8b)\begin{gather}\mathrm{i}\omega_0 E + \int v g \,{\rm d} v = j^{(S)}, \end{gather}

where $H = F/T$ and the variables $g$ and $E$ are now considered as the corresponding time-independent complex amplitudes. For definiteness, we impose a localized current source:

(2.9)\begin{equation} j^{(S)} = \mathrm{i}\omega_0 \exp({-(x-x_0)^2/(2 \varDelta_S^2)}). \end{equation}

Then, by (2.2), the corresponding charge density is that of an oscillating dipole:

(2.10)\begin{equation} \rho^{(S)} ={-}\frac{x-x_0}{\varDelta_S^2} \exp({-(x-x_0)^2/(2 \varDelta_S^2)}). \end{equation}

We will be interested in the spatial distribution of the electric field $E(x)$ driven by the source $j^{(S)}$ and undergoing linear Landau damping caused by interaction of this field with the distribution perturbation $g$.

To avoid numerical artifacts and keep the grid resolution reasonably low, we impose an artificial diffusivity $\eta$ in the velocity space by modifying (2.8a) as

(2.11)\begin{equation} \mathrm{i}\omega_0 g - v\partial_x g +\eta \partial_v^2g - v H E = 0. \end{equation}

This allows us to reduce the grid resolution in both velocity and real space.

To avoid numerical errors caused by the waves reflected from spatial boundaries, we impose outgoing (non-reflecting) boundary conditions on both edges (see Thompson Reference Thompson1987). The resulting system is

(2.12a)\begin{gather} \mathrm{i}\omega_0 g - \zeta^{\rm bc} v\partial_x g +\eta \partial_v^2g - v H E = 0, \end{gather}

(2.12b)\begin{gather}\mathrm{i}\omega_0 E + \int v g \,{\rm d} v = j^{(S)}, \end{gather}

where

(2.13)\begin{equation} \zeta^{\rm bc} = \left\{\begin{array}{@{}ll} 0 & \text{for incoming waves},\\ 1 & \text{otherwise}. \end{array}\right. \end{equation}

In our case of a one-dimensional system, the incoming waves correspond to $v < 0$ at the right spatial boundary of the simulation box and $v > 0$ at the left spatial boundary.

4.1 The vector ${\psi }$

First of all, to construct (1.1), one needs to encode $g_{j,k}$ and $E_j$ into ${\boldsymbol{\psi} }$. Since one needs to store $N_f = 2$ fields on a $N_x\times N_v$ phase space, the size of ${\boldsymbol{\psi} }$ should be $N_{\rm tot} = N_f N_{xv}$, where $N_{xv} = N_x N_v$, and this vector can be saved by using $1+n_x+n_v$ qubits, where $n_x$ and $n_v$ have been introduced in (3.3a,b). Within the vector ${\boldsymbol{\psi} }$, the fields are arranged in the following way:

(4.1)\begin{equation} \psi_{d N_{xv} + j N_v + k} = \left\{\begin{array}{@{}ll} g_{j,k}, & d = 0,\\ \delta_{k,0} E_j, & d = 1,\\ 0, & \text{otherwise}. \end{array}\right. \end{equation}

Here $\psi _{d N_{xv} + j N_v + k}$ are the elements of ${\boldsymbol{\psi} }$ with $j = [0, N_x)$, $k = [0,N_v)$ and $d = [0, N_f)$. Since the electric field does not depend on velocity, the second half of ${\boldsymbol{\psi} }$ is filled with $E_j$ only if the velocity index $k$ is equal to zero.

To address each field, $g_{j,k}$ or $E_j$, we introduce the register $r_f$ with one qubit. The zero state $\mid {0}\rangle_{r_f}$ corresponds to addressing the plasma distribution function, and the unit state $\mid {1}\rangle_{r_f}$ flags the electric field. We also use additional two registers, denoted $r_x$ and $r_v$, with $n_x$ and $n_v$ qubits, respectively, to specify the fields’ position in the real and velocity spaces, correspondingly. Then, one can express the vector ${\boldsymbol{\psi} }$ as

(4.2)\begin{equation} \mid{\psi}\rangle =\eta_{\psi, {\rm norm}} \sum_{j=0}^{N_x - 1} \left(\sum_{k=0}^{N_v-1} g_{j,k}\mid{0}\rangle_{r_f}\mid{j}\rangle_{r_x}\mid{k}\rangle_{r_v} + E_j\mid{1}\rangle_{r_f}\mid{j}\rangle_{r_x}\mid{0}\rangle_{r_v}\right), \end{equation}

where $\eta _{\psi, {\rm norm}}$ is the normalization factor used to ensure that $\langle \psi |\psi \rangle = 1$. (The assumed notation is such that the least significant qubit is the rightmost qubit. In quantum circuits, the least significant qubit is the lowest one.) Notably, this encoding can be called half-analogue, since the electric field and the electron distribution function are encoded into the amplitudes of the quantum state, which is a continuous variable. However, positions in phase space are discretized and digital since they are encoded into bitstrings of the quantum states.

4.2 The source term

Corresponding to (4.1), the discretized version of the vector $\boldsymbol {b}$ in (1.1) is as follows:

(4.3)\begin{equation} b_{d N_{xv} + j N_v + k} =\left\{\begin{array}{@{}ll} \delta_{k,0}j^{(S)}_j, & d = 1,\\ 0, & d = 0. \end{array}\right. \end{equation}

Note that all elements in the first half of $\boldsymbol {b}$ are zero, and, in its second half, only every $N_v$th element is non-zero. In the ‘ket’ notation, this can be written as

(4.4)\begin{equation} \mid{b}\rangle = \eta_{b,{\rm norm}}\sum_{j = 0}^{N_x-1} j^{(S)}_j\mid{1}\rangle_{r_f}\mid{j}\rangle_{r_x}\mid{0}\rangle_{r_v}, \end{equation}

where $\eta _{b,{\rm norm}}$ is the normalization factor used to ensure that $\langle b|b \rangle = 1$.

4.3 The matrix ${A}$

The corresponding $N_{\rm tot}\times N_{\rm tot}$ matrix $\boldsymbol {A}$ is represented as

(4.5)\begin{equation} \boldsymbol{A} = \begin{pmatrix} \boldsymbol{F} & \boldsymbol{C}^E \\ \boldsymbol{C}^f & \boldsymbol{S} \end{pmatrix}, \end{equation}

where $\boldsymbol {F}$, $\boldsymbol {C}^E$, $\boldsymbol {C}^f$ and $\boldsymbol {S}$ are $N_{xv}\times N_{xv}$ submatrices. A schematic structure of this matrix is shown in figure 1. The submatrix $\boldsymbol {F}$ encodes the coefficients in front of $g$ in the Vlasov equation (3.7) and is given by

(4.6)\begin{equation} \boldsymbol{F} = \begin{pmatrix} \boldsymbol{F}^{L,0} & \boldsymbol{F}^{L,1} & \boldsymbol{F}^{L,2} & 0 & 0 & \cdots \\ \boldsymbol{F}^{B,1} & \boldsymbol{F}^{B,0} & -\boldsymbol{F}^{B,1} & 0 & 0 & \cdots \\ 0 & \boldsymbol{F}^{B,1} & \boldsymbol{F}^{B,0} & -\boldsymbol{F}^{B,1} & 0 & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & 0 & \boldsymbol{F}^{B,1} & \boldsymbol{F}^{B,0} & -\boldsymbol{F}^{B,1} & 0 \\ \cdots & 0 & 0 & \boldsymbol{F}^{B,1} & \boldsymbol{F}^{B,0} & -\boldsymbol{F}^{B,1} \\ \cdots & 0 & 0 & \boldsymbol{F}^{R,2} & \boldsymbol{F}^{R,1} & \boldsymbol{F}^{R,0} \end{pmatrix}. \end{equation}

This submatrix consists of $N_x^2$ blocks, whose elements are mostly zeros. Each block is of size $N_v\times N_v$, and the position of each block's elements is determined by the row index $k_r = [0, N_v)$ and the column index $k_c = [0, N_v)$. The discretization at the left spatial boundary is described by the blocks $\boldsymbol {F}^{L,0}$, $\boldsymbol {F}^{L,1}$ and $\boldsymbol {F}^{L,2}$. The discretization at the right spatial boundary is described by the blocks $\boldsymbol {F}^{R,0}$, $\boldsymbol {F}^{R,1}$ and $\boldsymbol {F}^{R,2}$. The discretization at the bulk spatial points is described by the blocks $\boldsymbol {F}^{B,0}$ and $\boldsymbol {F}^{B,1}$. Using (3.8), (3.11) and (3.12), one obtains the following expressions for the elements in each block in (4.6) at the left spatial boundary:

(4.7a)\begin{gather} F^{L,0}_{k_r, k_c} = \delta_{k_r, k_c} P_{0,k_r} + P^{(v)}_{k_c, k_r}, \end{gather}

(4.7b)\begin{gather}F_{k_r, k_c}^{L,1} ={-}4 v_{k_r} \sigma \zeta^{\rm bc}_{0,k_r} \delta_{k_r, k_c}, \end{gather}

(4.7c)\begin{gather}F_{k_r, k_c}^{L,2} = v_{k_r} \sigma \zeta^{\rm bc}_{0,k_r} \delta_{k_r, k_c}. \end{gather}

At the right spatial boundary,

(4.8a)\begin{gather} F^{R,0}_{k_r, k_c} = \delta_{k_r, k_c} P_{q_x,k_r} + P^{(v)}_{k_c, k_r}, \end{gather}

(4.8b)\begin{gather}F_{k_r, k_c}^{R,1} = 4 v_{k_r} \sigma \zeta^{\rm bc}_{q_x,k_r} \delta_{k_r, k_c}, \end{gather}

(4.8c)\begin{gather}F_{k_r, k_c}^{R,2} ={-}v_{k_r} \sigma \zeta^{\rm bc}_{q_x,k_r} \delta_{k_r, k_c}, \end{gather}

and at bulk spatial points,

(4.9a)\begin{gather} F^{B,0}_{k_r, k_c} = \delta_{k_r, k_c} P_{j,k_r} + P^{(v)}_{k_c, k_r}, \end{gather}

(4.9b)\begin{gather}F_{k_r, k_c}^{B,1} = v_{k_r} \sigma \delta_{k_r, k_c}. \end{gather}

We denote the part of the submatrix $\boldsymbol {F}$ that depends on $v$ as $\tilde {\boldsymbol {F}}$. The matrix elements of $\tilde {\boldsymbol {F}}$ are indicated in figure 1 with red. The part of $\boldsymbol {F}$ that does not depend on $v$ is denoted as $\hat {\boldsymbol {F}}$ and shown in blue in figure 1 within the submatrix $\boldsymbol {F}$.

Figure 1. A schematic showing the structure of the matrix $\boldsymbol {A}$ (4.5) with $n_x = 3$ and $n_v = 3$. The solid lines separate the submatrices introduced in (4.5). The dashed lines indicate blocks of size $N_v\times N_v$. The blue markers indicate velocity-independent elements.

The matrix $\boldsymbol {C}^E$ is block diagonal and encodes the coefficients in front of $E$ in the Vlasov equation (3.7):

(4.10)\begin{equation} C^E_{j_r N_v + k_r, j_c N_v + k_c} ={-}\delta_{j_r,j_c}\delta_{k_c,0}v_{k_r} H_{j_r,k_r}. \end{equation}

Here $j_r,j_c = [0, N_x)$ and $k_r, k_c = [0,N_v)$. In this submatrix the first column in each diagonal block of size $N_v\times N_v$ is non-sparse whilst all other columns are filled with zeros. The Kronecker delta $\delta _{k_c,0}$ appears because of the chosen encoding of the electric field into the state vector ${\boldsymbol{\psi} }$ according to (4.1).

The submatrix $\boldsymbol {C}^f$ is also block diagonal and encodes the coefficients in front of $g$ in (3.13):

(4.11)\begin{equation} C^f_{j_r N_v + k_r, j_c N_v + k_c} = \delta_{j_r,j_c}\delta_{k_r,0}v_{k_c}. \end{equation}

Here $j_r,j_c = [0, N_x)$ and $k_r, k_c = [0,N_v)$. The first row in each block of size $N_v\times N_v$ is non-sparse due to the sum in (3.13).

Finally, the matrix $\boldsymbol {S}$ is diagonal and encodes the coefficients in front of $E$ in (3.13):

(4.12)\begin{equation} S_{j_r N_v + k_r, j_c N_v + k_c} = \delta_{j_r,j_c}\delta_{k_r,k_c}\mathrm{i}\omega_0. \end{equation}

Here $j_r,j_c = [0, N_x)$ and $k_r, k_c = [0,N_v)$.

6.1 Initialization

To encode the right-hand-side vector (4.4) into a quantum circuit, one can use the fact that the shape of the source current (2.9) is Gaussian. As shown in Novikau, Dodin & Startsev (Reference Novikau, Dodin and Startsev2023), Kane, Gomes & Kreshchuk (Reference Kane, Gomes and Kreshchuk2023) and Hariprakash et al. (Reference Hariprakash, Modi, Kreshchuk, Kane and Bauer2023), one can encode this function by using either QSVT (see Gilyén et al. Reference Gilyén, Su, Low and Wiebe2019; Martyn et al. Reference Martyn, Rossi, Tan and Chuang2021) or the so-called quantum eigenvalue transformation of unitaries (QETU) (see Dong, Lin & Tong Reference Dong, Lin and Tong2022), where the scaling of the resulting circuit is $\mathcal {O}(n_v \log _2(\varepsilon ^{-1}_{\rm qsvt}))$ and $\varepsilon _{\rm qsvt}$ is the desired absolute error in the QSVT approximation of the Gaussian. However, one should keep in mind that the success probability of the initialization circuit depends on the Gaussian width. To increase this probability, amplitude amplification can be used (see Brassard et al. Reference Brassard, Høyer, Mosca and Tapp2002). After the initialization, the required initial state $\mid {b}\rangle$ is usually entangled with the zero state of the ancillae used for the QSVT/QETU procedure and for the amplitude amplification. Thus, the subsequent QSVT circuit computing the inverse matrix $\boldsymbol {A}^{-1}$ should be controlled by this zero state to guarantee that $\boldsymbol {A}^{-1}$ acts on the state $\mid {b}\rangle$.

6.2 Block encoding: basic idea

To solve (1.1) with the matrix (4.5) and the source (4.4) on a digital quantum computer, one can use the QSVT, which approximates the inverse matrix $\boldsymbol {A}^{-1}$ with an odd polynomial of the matrix singular values. (The mathematical foundations of the QSVT are described, for example, in Gilyén et al. Reference Gilyén, Su, Low and Wiebe2019; Martyn et al. Reference Martyn, Rossi, Tan and Chuang2021; Lin Reference Lin2022.) The QSVT returns a quantum state $\mid {\psi _{\rm qsvt}}\rangle$ whose projection on zero ancillae is proportional to the solution ${\boldsymbol{\psi} }$ of (1.1):

(6.1)\begin{equation} \mid{\psi_{\rm qsvt}}\rangle = \frac{\exp({\mathrm{i}\phi_{\rm glob}})}{\kappa_{\rm qsvt}}\boldsymbol{A}^{{-}1} \mid{0}\rangle_{\rm anc}\mid{b}\rangle_{\rm in} + \mid{\neq 0}\rangle_{\rm anc}\mid{\ldots}\rangle_{\rm in}. \end{equation}

Here $\phi _{\rm glob}$ is an unknown global angle and the scalar parameter $\kappa _{\rm qsvt}$ is of the order of the condition number of the matrix $\boldsymbol {A}$.

A typical QSVT circuit is shown in figure 7, where the angles $\phi _i$ are computed classically. These angles serve as the parameters that specify the function computed by the QSVT circuit. In our case, the function is the inverse function. (More details about the computation of the QSVT angles can be found in Dong et al. Reference Dong, Meng, Whaley and Lin2021; Ying Reference Ying2022). The subcircuit $U_A$ is the so-called BE oracle with the following matrix representation:

(6.2)\begin{equation} \boldsymbol{U}_A =\begin{pmatrix} \boldsymbol{A} & \cdot \\ \cdot & \cdot \end{pmatrix}. \end{equation}

Here $\boldsymbol {U}_A$ is a unitary matrix. (The dots correspond to submatrices that keep $\boldsymbol {U}_A$ unitary but otherwise are unimportant.) This unitary encodes the matrix $\boldsymbol {A}$ as a sub-block that is accessed by setting the ancilla register ‘be’ to zero:

(6.3)\begin{equation} \boldsymbol{U}_A\mid{0}\rangle_{\rm be}\mid{\psi}\rangle_{\rm in} = \boldsymbol{A}\mid{0}\rangle_{\rm be} \mid{\psi}\rangle_{\rm in} + \mid{\neq 0}\rangle_{\rm be}\mid{\ldots}\rangle_{\rm in}. \end{equation}

The QSVT addresses the BE oracle $\mathcal {O}(\kappa _{\rm qsvt}\log _2(\varepsilon ^{-1}_{\rm qsvt}))$ times to approximate the inverse matrix $\boldsymbol {A}^{-1}$. The efficient implementation of the BE oracle is the key to a potential quantum speedup that might be provided by the QSVT. As discussed in Novikau et al. (Reference Novikau, Dodin and Startsev2023), the QSVT can provide a polynomial speedup for two- and higher-dimensional classical wave systems if the quantum circuit of the BE oracle scales polylogarithmically or better with the size of $\boldsymbol {A}$.

Figure 7. The QSVT circuit encoding a real polynomial of order $N_a$, where $N_a$ is odd, by using $N_a + 1$ angles $\phi _k$ pre-computed classically. The gates denoted as $R_{z,k}$ represent the rotations $R_z(2\phi _k)$. For an even $N_a$, the gate $Z$ should be removed and the rightmost BE oracle $U_A$ should be replaced with its Hermitian adjoint version $U_A^{\dagger}$.

To block encode a non-Hermitian matrix such as $\boldsymbol {A}$, one can first extend it to a Hermitian one as

(6.4)\begin{equation} \boldsymbol{A}_{\rm ext} = \begin{pmatrix} 0 & \boldsymbol{A} \\ \boldsymbol{A}^{\dagger} & 0 \end{pmatrix} \end{equation}

and then use the technique from Berry & Childs (Reference Berry and Childs2012). However, this will require at least two additional ancilla qubits for the extension (6.4). Another option is to decompose $\boldsymbol {A}$ into two Hermitian matrices, i.e.

(6.5)\begin{equation} \boldsymbol{A} = \boldsymbol{A}_h + \mathrm{i} \boldsymbol{A}_a, \end{equation}

where $\boldsymbol {A}_h = (\boldsymbol {A} + \boldsymbol {A}^{\dagger} )/2$ and $\boldsymbol {A}_a = (\boldsymbol {A} - \boldsymbol {A}^{\dagger} )/(2\mathrm {i})$. The sum (6.5) can be computed by using the circuit of linear combination of unitaries, which requires an additional ancilla. Also note that the matrices $\boldsymbol {A}_a$ and $\boldsymbol {A}_h$, although being Hermitian, still have a non-trivial structure. Thus, for this method, it is necessary to block encode two separate matrices, which may double the depth of the BE oracle.

To reduce the number of ancillae and avoid encoding two matrices instead of one, we propose to encode the non-Hermitian $\boldsymbol {A}$ directly, without invoking the extension (6.4) or the splitting (6.5). This technique was already used in Novikau et al. (Reference Novikau, Dodin and Startsev2023). Although the direct encoding requires ad hoc construction of some parts of the BE oracle, this approach leads to a more compact BE circuit.

To encode $\boldsymbol {A}$ into the unitary $\boldsymbol {U}_A$, one needs to normalize $\boldsymbol {A}$ such that $\varsigma ||\boldsymbol {A}||_{\rm max} \leq 1$, where

(6.6)\begin{equation} ||\boldsymbol{A}||_{\rm max} = \max_k\sum_j\sqrt{|A_{kj}|^2} \end{equation}

and $\varsigma$ is related to the matrix non-sparsity as will be explained later (see (6.14)). Here, by the term ‘non-sparsity’, we understand the maximum number of non-zero elements in a matrix row or column.Footnote ³ Hence, $\boldsymbol {A}$ should be normalized as

(6.7)\begin{equation} \boldsymbol{A} \to \boldsymbol{A}/\left(||\boldsymbol{A}||_{\rm max}\varsigma\right). \end{equation}

The general structure of a BE oracle implementing the direct encoding of a non-Hermitian matrix is

(6.8)\begin{equation} U_A = O^{\dagger}_F O^{B}_{F,{\rm corr}} O_M O_H O^{F}_{F,{\rm corr}} O_F, \end{equation}

where the oracles $O_F$, $O^{F}_{F,{\rm corr}}$, $O_M$, $O^{B}_{F,{\rm corr}}$ and $O^{\dagger} _F$ encode the positions of non-zero matrix elements in $\boldsymbol {A}$, and $O_H$ encodes the values of these elements. The upper superscripts ‘F’ and ‘B’ stand for ‘forward’ and ‘backward’ action. The oracle $U_A$ can be constructed by using quantum gates acting on a single qubit but controlled by multiple qubits. Here, such gates are called single-target multicontrolled (STMC) gates (§ 6.3).

Equation (6.8) is based on the BE technique from Berry & Childs (Reference Berry and Childs2012) with the only difference that the oracles $O^{F}_{F,{\rm corr}}$ and $O^{B}_{F,{\rm corr}}$ are introduced to take into account the non-Hermiticity of $\boldsymbol {A}$ and are computed ad hoc by varying STMC gates. Basically, the oracles $O_F$, $O_M$ and $O^{\dagger} _F$ create a structure of a preliminary Hermitian matrix that is as close as possible to the target non-Hermitian matrix $\boldsymbol {A}$. The structure is then corrected by the oracles $O^{F}_{F,{\rm corr}}$ and $O^{B}_{F,{\rm corr}}$ constructed by varying their circuits to encode the structure of $\boldsymbol {A}$.

We consider separately the action of the oracle

(6.9a)\begin{align} U_D = O^{B}_F O_M O^{F}_F, \end{align}

(6.9b)\begin{align} O^{F}_F = O^{F}_{F,{\rm corr}} O_F, \end{align}

(6.9c)\begin{align} O^{B}_F = O^{\dagger}_F O^{B}_{F,{\rm corr}}, \end{align}

where the oracle $O_H$ is not included, and the oracles $O^{F}_F$ and $O^{B}_F$ are introduced to simplify notations. The matrix representation of the oracle $U_D$ projected onto zero ancillae is denoted as $\boldsymbol {D}^A$. The matrix $\boldsymbol {D}^A$ has non-zero elements at the same positions as those of the non-zero elements of $\boldsymbol {A}$.

We use the input register ‘in’ to encode a row index $i_r$ and the ancilla register $a_c$ to perform intermediate computations. Then,

(6.10)\begin{equation} U_D\mid{0}\rangle_{a_c}\mid{i_r}\rangle_{\rm in} = \sum_{i_c = 0}^{N_{c,i_r}-1} d_{i_r i_{c}} \mid{0}\rangle_{a_c}\mid{i_{c}}\rangle_{\rm in} + \mid{\neq 0}\rangle_{a_c}\mid{\ldots}\rangle_{\rm in}, \end{equation}

where $i_{c}$ are the column indices of all non-zero elements of $\boldsymbol {A}$ at the row $i_r$; $N_{c,i_r}$ is the number of these elements at $i_r$; $d_{i_r i_{c}} \leq 1$ are the matrix elements of $\boldsymbol {D}^A$. The number $N_{c,i_r}$ is less or equal to the non-sparsity of $\boldsymbol {A}$ and can be different at different $i_c$. The elements $d_{i_r i_{c}}$ are usually different powers of the factor $2^{-1/2}$ that appears due to the usage of multiple Hadamard gates $H$.

The oracle $O_H$, which enters $U_A$ but is not a part of $U_D$, takes the matrix $\boldsymbol {D}^A$ and modifies its elements to form $\boldsymbol {A}$. Usually, $O_H$ acts on an extra ancilla $a_e$ that is not used by $U_D$. Due to that, we can formally consider $U_D$ separately from $O_H$. The action of $O_H$ can be considered as the mapping

(6.11)\begin{equation} O_H: d_{i_r i_c} \to A_{i_r i_c}, \end{equation}

where $A_{i_r i_c}$ are elements of the matrix $\boldsymbol {A}$ after the normalization (6.7). For instance, to encode a real-value element $A_{i_r i_c}$, one can use the rotation gate $R_y(\theta )$:

(6.12)\begin{equation} R_y(\theta) d_{i_r i_c}\mid{0}\rangle_{a_e} = \cos(\theta/2)d_{i_r i_c}\mid{0}\rangle_{a_e} + \ldots\mid{1}\rangle_{a_e}. \end{equation}

The factor $d_{i_r i_c}$ appears from the action of the oracle $U_D$ (6.10). Our goal is to have $A_{i_r i_c} = \cos (\theta /2)d_{i_r i_c}$. Thus,

(6.13)\begin{equation} \theta = 2 \arccos(A_{i_r i_c} / d_{i_r i_c}). \end{equation}

The fact that $d_{i_r i_c} \leq 1$ is the reason why it is necessary to include $\varsigma$ into the normalization (6.7). From this, we conclude that

(6.14)\begin{equation} \varsigma = \max_{i,j, D^A_{ij}\neq 0}|(D^A_{ij})^{{-}1}|. \end{equation}

The oracle $O_H$ usually consists of STMC rotations gates $R_x$, $R_y$, $R_z$ and $R_c$ (C1). The first two are used to encode imaginary and real values, correspondingly. The third one can be used to change the sign of a value or to turn a real value into an imaginary one if necessary, and vice versa. The gate $R_c$ is used to encode complex values.

Now, let us specify the action of different parts of $U_A$. The oracles $O^{F}_{F,{\rm corr}} O_F$ encode column indices into the ancilla register $a_c$:

(6.15)\begin{equation} O^{F}_F\mid{0}\rangle_{a_e}\mid{0}\rangle_{a_c}\mid{i_r}\rangle_{\rm in} = \sum_{i_c}\sqrt{d_{i_r i_c}}\mid{0}\rangle_{a_e}\mid{i_c}\rangle_{a_c}\mid{i_r}\rangle_{\rm in}. \end{equation}

The oracle $O_H$ uses the row index from the state register ‘in’ and the column indices from the ancilla register $a_c$ to determine which element should be computed and then encodes it into the state amplitude:

(6.16)\begin{equation} O_H \sum_{i_c}\sqrt{d_{i_r i_c}}\mid{0}\rangle_{a_e}\mid{i_c}\rangle_{a_c}\mid{i_r}\rangle_{\rm in} = \sum_{i_c}\frac{A_{i_r i_c}}{\sqrt{d_{i_r i_c}}}\mid{0}\rangle_{a_e}\mid{i_c}\rangle_{a_c}\mid{i_r}\rangle_{\rm in} + \mid{\neq 0}\rangle_{a_e}\mid{\ldots}\rangle. \end{equation}

After that, the oracle $O_M$ transfers the column indices from $a_c$ to the input register:

(6.17)\begin{align} & O_M\left(\sum_{i_c}\frac{A_{i_r i_c}}{\sqrt{d_{i_r i_c}}}\mid{0}\rangle_{a_e}\mid{i_c}\rangle_{a_c}\mid{i_r}\rangle_{\rm in} + \mid{\neq 0}\rangle_{a_e}\mid{\ldots}\rangle\right) \nonumber\\ & \quad = \sum_{i_c}\frac{A_{i_r i_c}}{\sqrt{d_{i_r i_c}}}\mid{0}\rangle_{a_e}\mid{i_r}\rangle_{a_c}\mid{i_c}\rangle_{\rm in} + \mid{\neq 0}\rangle_{a_e}\mid{\ldots}\rangle. \end{align}

Finally, the oracles $O^{\dagger} _F O^{B}_{F,{\rm corr}}$ entangle the states encoding the column indices in the input register with the zero state in the ancilla register:

(6.18)\begin{align} & O^{B}_F\left(\sum_{i_c}\frac{A_{i_r i_c}}{\sqrt{d_{i_r i_c}}}\mid{0}\rangle_{a_e} \mid{i_r}\rangle_{a_c}\mid{i_c}\rangle_{\rm in}+ \mid{\neq 0}\rangle_{a_e}\mid{\ldots}\rangle\right) \nonumber\\ & \quad = \sum_{i_c}A_{i_r i_c}\mid{0}\rangle_{a_e}\mid{0}\rangle_{a_c}\mid{i_c}\rangle_{\rm in} + \mid{\neq 0}\rangle_{a_e}\mid{\neq 0}\rangle_{a_c}\mid{\ldots}\rangle. \end{align}

6.3 Single-target multicontrolled gates

If one has a single-target gate $G$, whose matrix representation is

(6.19)\begin{equation} \boldsymbol{G} = \begin{pmatrix} G_{00} & G_{01} \\ G_{10} & G_{11} \end{pmatrix}, \end{equation}

then the corresponding STMC gate $C_{\{q_{c\delta }\}}G^{(q_t)}$ is defined as the gate $G$ acting on the target qubit $q_t$ and controlled by a set of qubits $\{q_{c\delta }\}$. If $q_{c\delta }$ is a control qubit then the gate $G$ is triggered if and only if $\mid {\delta }\rangle_{q_{c\delta }}$, where $\delta = 0$ or $1$. If the gate $G$ acting on a quantum state vector ${\boldsymbol{\psi} }$ of $n$ qubits is controlled by the qubit $q_{c\delta } \in [0,n)$, then only the state vector's elements with the indices $\{i_e\}_{q_{c\delta }}$ can be modified by the gate $G$:

(6.20a)\begin{gather} i_e = 2 N_{\rm ctrl} j_b + j_{\rm step} + \delta N_{\rm ctrl}, \end{gather}

(6.20b)\begin{gather}N_{\rm ctrl} = 2^{q_{c\delta}}, \end{gather}

(6.20c)\begin{gather}j_b = [0,-1+2^n/(2N_{\rm ctrl})], \end{gather}

(6.20d)\begin{gather}j_{\rm step} = [0,N_{\rm ctrl}). \end{gather}

For instance, for $q_{c1} = 0$, every second element of ${\boldsymbol{\psi} }$ can be modified by the gate $G$. Then, the STMC gate $C_{\{q_{c\delta }\}}G^{(q_t)}$ can modify only the elements from the set

(6.21)\begin{equation} \mathcal{S} = \bigcap_{q_{c\delta,k}\in \{q_{c\delta}\}} \{i_e\}_{q_{c\delta,k}}, \end{equation}

where $\bigcap$ is the intersection operator. The most common case is when $q_t$ is a more significant qubit than the control ones, and the initial state of $q_t$ is the zero state. In this case, the action of $C_{\{q_{c\delta }\}}G^{(q_t)}$ can be described as

(6.22)\begin{align} & C_{\{q_{c\delta}\}}G^{(q_t)}\mid{0}\rangle_{q_t} \sum_{k\in\{i_c\}_{\rm init}} \eta_{k}\mid{k}\rangle_{\rm ctrl} \nonumber\\ & \quad = (G_{00}\mid{0}\rangle_{q_t} + G_{10}\mid{1}\rangle_{q_t}) \sum_{k\in \{i_c\}_{\rm init}\cap \mathcal{S}} \eta_{k}\mid{k}\rangle_{\rm ctrl} + \mid{0}\rangle_{q_t}\sum_{k\in \{i_c\}_{\rm init}\setminus\mathcal{S}} \eta_{k}\mid{k}\rangle_{\rm ctrl}, \end{align}

where $\eta _{k}$ are the complex amplitudes of the initial state of the control register.

6.4 General algorithm for the BE

To construct the BE oracle $U_A$, we use the following general procedure.

1. (i) Introduce ancilla qubits necessary for intermediate computations in the oracle $U_A$ (§ 6.5).

2. (ii) Assume that the bitstring of the qubits of the input register (also called the state register) encodes a row index of the matrix $\boldsymbol {A}$.

3. (iii) Using STMC gates, construct the oracle $U_D$ following the idea presented in (6.10) to encode column indices as a superposition of bitstrings in the state register (§ 6.6).

4. (iv) Compute the matrix $\boldsymbol {D}^A$ or derive it using several matrices $\boldsymbol {D}^A$ constructed for matrices $\boldsymbol {A}$ of small sizes (§ 7.3).

5. (v) Normalize the matrix $\boldsymbol {A}$ according to (6.7) using the non-sparsity-related parameter (6.14).

6. (vi) Using STMC rotation gates, construct the oracle $O_H$ to perform the transformation (6.11) (§ 7).

Once the circuit for the oracle $U_A$ is constructed using STMC gates, one can transpile the circuit into a chosen universal set of elementary gates. Standard decomposition methods require at least $\mathcal {O}(n)$ of basic gates (see Barenco et al. Reference Barenco, Bennett, Cleve, DiVincenzo, Margolus, Shor, Sleator, Smolin and Weinfurter1995), where $n$ is the number of controlling qubits in an STMC gate. Yet, it was recently shown (see Claudon et al. Reference Claudon, Zylberman, Feniou, Debbasch, Peruzzo and Piquemal2023) that it is possible to decompose an arbitrary STMC gate into a circuit with $\mathcal {O}(\log _2(n)^{\log _2(12)}\log _2(1/\epsilon _{\rm STMC}))$ depth where $\epsilon _{\rm STMC}$ is the allowed absolute error in the approximation of the STMC gate. In our assessment of the BE oracle's scaling below, we assume that the corresponding circuit comprises STMC gates not decomposed into elementary gates.

6.5 Ancilla qubits for the BE

The first step to construct a circuit of the BE oracle is to introduce ancilla qubits and assign the meaning to their bitstrings. As seen from figure 1, the matrix $\boldsymbol {A}$ is divided into four submatrices. To address each submatrix, we introduce the ancilla qubit $a_f$. In combination with the register $r_f$ introduced in (4.2), the register $a_f$ allows addressing each submatrix in the following way:

(6.23a)\begin{gather} \mid{0}\rangle_{a_f}\mid{0}\rangle_{r_f}\to \boldsymbol{F}, \end{gather}

(6.23b)\begin{gather}\mid{1}\rangle_{a_f}\mid{0}\rangle_{r_f}\to \boldsymbol{C}^E, \end{gather}

(6.23c)\begin{gather}\mid{0}\rangle_{a_f}\mid{1}\rangle_{r_f}\to \boldsymbol{C}^f, \end{gather}

(6.23d)\begin{gather}\mid{1}\rangle_{a_f}\mid{1}\rangle_{r_f}\to \boldsymbol{S}. \end{gather}

If $i_r$ and $i_c$ are row and column indices of $\boldsymbol {A}$, respectively, such that

(6.24a)\begin{gather} i_r = i_{fr} N_{xv} + i_{xr}N_v + i_{vr}, \end{gather}

(6.24b)\begin{gather}i_c = i_{fc} N_{xv} + i_{xc}N_v + i_{vc}, \end{gather}

with integers $i_{fr}, i_{fc} = [0,1]$, $i_{xr}, i_{xc} = [0, N_x)$ and $i_{vr}, i_{vc} = [0, N_v)$, then the ancilla $a_f$ encodes the integer $i_{fc}$.

According to figure 1, each submatrix is split into $N_x^2$ blocks of size $N_v\times N_v$ each. Within each submatrix, only the diagonal and a few off-diagonal (in the submatrix $\boldsymbol {F}$) blocks contain non-zero matrix elements. To address these blocks, we introduce the ancilla register $a_{xr}$ with three qubits. This register stores the relative positions of blocks with respect to the main diagonal within each submatrix. The states $\mid {000}\rangle_{a_{xr}}$ and $\mid {100}\rangle_{a_{xr}}$ correspond to the blocks on the main diagonal of a submatrix. The states $\mid {001}\rangle_{a_{xr}}$ and $\mid {010}\rangle_{a_{xr}}$ correspond to the off-diagonal blocks shifted by one and two blocks to the right, respectively. The states $\mid {101}\rangle_{a_{xr}}$ and $\mid {110}\rangle_{a_{xr}}$ correspond to the off-diagonal blocks shifted by one and two blocks to the left, respectively. Using the same notation as in (6.24), the meaning of the register $a_{xr}$ can be described schematically as follows:

(6.25a)\begin{gather} \mid{000}\rangle_{a_{xr}}\ \text{and}\ \mid{100}\rangle_{a_{xr}} \to i_{xc} = i_{xr}, \end{gather}

(6.25b)\begin{gather}\mid{001}\rangle_{a_{xr}} \to i_{xc} = i_{xr} + 1, \end{gather}

(6.25c)\begin{gather}\mid{010}\rangle_{a_{xr}} \to i_{xc} = i_{xr} + 2, \end{gather}

(6.25d)\begin{gather}\mid{101}\rangle_{a_{xr}} \to i_{xc} = i_{xr} - 1, \end{gather}

(6.25e)\begin{gather}\mid{110}\rangle_{a_{xr}} \to i_{xc} = i_{xr} - 2. \end{gather}

To encode the above bitstrings, one can use the circuits shown in figures 8(a) and 8(b).

Figure 8. (a) The circuit encoding the superposition of states $\mid {000}\rangle_q$, $\mid {001}\rangle_q$ and $\mid {101}\rangle_q$ for a given input zero state. (b) The circuit encoding the superposition of states $\mid {000}\rangle_q$, $\mid {001}\rangle_q$ and $\mid {010}\rangle_q$ for a given input state produced by the circuit 8(a). To encode the superposition of states $\mid {100}\rangle_q$, $\mid {101}\rangle_q$ and $\mid {110}\rangle_q$, one uses the same circuit except that an additional $X$ gate is applied to the qubit $q_2$ in the end. (c) The circuit encoding the superposition of states $\mid {000}\rangle_q$, $\mid {001}\rangle_q$, $\mid {010}\rangle_q$ and $\mid {011}\rangle_q$ for a given input state produced by the circuit 8(a). To encode the superposition of states $\mid {100}\rangle_q$, $\mid {101}\rangle_q$, $\mid {110}\rangle_q$ and $\mid {111}\rangle_q$, one uses the same circuit except that an additional $X$ gate is applied to the qubit $q_2$ in the end.

In the submatrix $\boldsymbol {C}^f$, some of the blocks have rows that contain $N_v$ non-zero elements. To address these elements, we introduce the ancilla register $a_v$ with $n_v$ qubits. This register stores the integer $i_{vc}$ (6.24) when the elements of $\boldsymbol {C}^f$ and $\boldsymbol {C}^E$ are addressed; otherwise, the register is not used. For instance, to address the non-zero elements in the submatrix $\boldsymbol {C}^f$, one uses the following encoding:

(6.26)\begin{equation} \mid{0}\rangle_{a_f}\mid{i_{vc}}\rangle_{a_v}\mid{1}\rangle_{r_f}\mid{i_x}\rangle_{r_x} \mid{0}\rangle_{r_v} \to C^f_{i_x N_v,\ i_x N_v + i_{vc}}. \end{equation}

As regards to $\boldsymbol {C}^E$, since all its non-zero elements have $i_{vc} = 0$, one keeps the register $a_v$ in the zero state when $\mid {1}\rangle_{a_f}\mid {0}\rangle_{r_f}$.

The ancilla register $a_{vr}$ with three qubits is introduced to encode positions of the matrix elements within the non-zero blocks of the submatrix $\boldsymbol {F}$. In particular, the states $\mid {000}\rangle_{a_{vr}}$ and $\mid {100}\rangle_{a_{vr}}$ correspond to the elements on the main diagonal of a block. The states $\mid {001}\rangle_{a_{vr}}$, $\mid {010}\rangle_{a_{vr}}$ and $\mid {011}\rangle_{a_{vr}}$ correspond to the elements shifted by one, two and three cells, respectively, to the right from the diagonal. The states $\mid {101}\rangle_{a_{vr}}$, $\mid {110}\rangle_{a_{vr}}$ and $\mid {111}\rangle_{a_{vr}}$ correspond to the elements shifted by one, two and three cells, respectively, to the left from the diagonal. Using the notations from (6.24), the meaning of the register $a_{vr}$ can be described as

(6.27a)\begin{gather} \mid{000}\rangle_{a_{vr}}\ \text{and}\ \mid{100}\rangle_{a_{vr}} \to i_{vc} = i_{vr}, \end{gather}

(6.27b)\begin{gather}\mid{001}\rangle_{a_{vr}} \to i_{vc} = i_{vr} + 1, \end{gather}

(6.27c)\begin{gather}\mid{010}\rangle_{a_{vr}} \to i_{vc} = i_{vr} + 2, \end{gather}

(6.27d)\begin{gather}\mid{011}\rangle_{a_{vr}} \to i_{vc} = i_{vr} + 3, \end{gather}

(6.27e)\begin{gather}\mid{101}\rangle_{a_{vr}} \to i_{vc} = i_{vr} - 1, \end{gather}

(6.27f)\begin{gather}\mid{110}\rangle_{a_{vr}} \to i_{vc} = i_{vr} - 2, \end{gather}

(6.27g)\begin{gather}\mid{111}\rangle_{a_{vr}} \to i_{vc} = i_{vr} - 3. \end{gather}

To encode the above bitstrings, one can use the circuits shown in figures 8(b) and 8(c).

Also, we introduce the ancilla $a_e$ whose zero-state's amplitude will encode the complex value of a given matrix element.

6.6 Constructing the oracle ${U_D}$

A schematic of the BE oracle $U_A$ is shown in figure 9, where the dashed blue box is the oracle $O_H$ described in § 7, and the rest of the boxes are components of the oracle $U_D$.

Figure 9. The circuit representation of the BE oracle $U_A$. The oracle $O_M$ is shown in figure 11. The oracles $F$, $C^f$, $C^E$ and $S$ encode the corresponding submatrices introduced in (4.5), and are described in § 7. Here, the symbol $\varnothing$ means that the gate does not use the corresponding qubit. The qubit $r_{\rm sin}$ shown in some oracles indicates that the corresponding oracle includes the circuit shown in figure 15. The qubit $r_{\rm qsvt}$ indicates that the oracle $C^E$ is constructed using QSVT. The positions of the qubits $r_f$ and $a_f$ are changed for the sake of clarity.

The oracle $U_D$ is built ad hoc manually by considering its suboracles separately. The suboracle $O_F^{F} = O^{F}_{F,{\rm corr}} O_F$ is constructed based on (6.15), where the ‘in’ register includes the registers $r_v$, $r_x$ and $r_f$ described in § 4.1. The ancilla registers $a_v$, $a_{vr}$, $a_{xr}$ and $a_f$ described in § 6.5 represent the $a_c$ register in (6.15). For instance, the part of the circuit $O_F$ encoding the ancilla $a_f$ (6.23) is shown in figure 10, where the first Hadamard gate creates the address to the submatrices $\boldsymbol {F}$ and $\boldsymbol {C}^E$, and the Pauli $X$ gate generates the address to the elements of the submatrix $\boldsymbol {S}$ with $i_{vr} > 0$ (using the notation from (6.24)). The last Hadamard gate produces the addresses to $\boldsymbol {C}_f$ and $\boldsymbol {S}$ at $i_{vr} = 0$.

Figure 10. The circuit to encode information into the ancilla $a_f$ according to (6.23). The registers $r_v$, $r_x$ and $r_f$ encode the row index $i_r$, i.e. the indices $i_{vr}$, $i_{xr}$ and $i_{fr}$, correspondingly, according to (6.24).

The implementation of the suboracle $O_F^{B} = O^{\dagger} _F O^{B}_{F,{\rm corr}}$ is done based on (6.18). Its correcting part $O^{B}_{F,{\rm corr}}$ is implemented ad hoc by varying positions and control nodes of STMC gates.

The suboracle $O_M$ that performs the mapping (6.17) uses the SWAP gates to exchange the states in the registers $a_f$ and $a_v$ with the states in the registers $r_f$ and $r_v$, respectively, as shown in figure 11. To encode absolute column indices into the input registers using the states of the registers $a_{xr}$ and $a_{vr}$, one should follow the rules (6.25) and (6.27) and apply the quantum adders and subtractors described in Appendix C. The important feature of the oracle $O_M$ is that the number of arithmetic operators in it does not depend on $N_x$ or $N_v$. Since the circuits of the arithmetic operators scale as $\mathcal {O}(n_x)$ or $\mathcal {O}(n_v)$ depending on the target register of the operator, the scaling of $O_M$ is $\mathcal {O}(n_x + n_v)$. The full circuit of $U_D$ can be found in Novikau (Reference Novikau2024b). The modelling of the circuit is done using the QuCF framework (see Novikau Reference Novikau2024c), which is based on the QuEST toolkit (see Jones et al. Reference Jones, Brown, Bush and Benjamin2019).

Figure 11. The circuit of the oracle $O_M$. The adders $A1$–$A3$ and the subtractors $S1$–$S3$ are described in Appendix C.

7.1 General strategy

The purpose of the oracle $O_H$ is to compute $A_{i_r i_c}$ following (6.13). Therefore, further in the text, we consider the following rescaled matrix elements:

(7.1)\begin{equation} A_{i_r i_c} \to \frac{A_{i_r i_c}}{|d_{i_r i_c}|}. \end{equation}

One can see from § 4.3 that two types of elements can be distinguished in $\boldsymbol {A}$. The first type includes elements that depend only on the discretization and system parameters such as the antenna frequency $\omega _0$ and the sizes of the grid cells. These elements sit in the submatrix $\boldsymbol {S}$ (4.12) and the submatrix $\hat {\boldsymbol {F}}$ introduced after (4.9). Since these elements appear mainly from the spatial, velocity and time derivatives, they remain mostly constant at bulk spatial and bulk velocity points, but become highly intermittent at spatial and velocity boundaries. The values of the elements of $\hat {\boldsymbol {F}}$ strongly vary at boundaries of each $N_v\times N_v$ block. Note that the submatrix $\hat {\boldsymbol {F}}$, which is a part of the matrix $\boldsymbol {A}$, is also subject to the rescaling (7.1). Thus, its elements change from one block to another depending not only on the original values of $\hat {\boldsymbol {F}}$ (4.7), (4.8) and (4.9) but also on the specific implementation of $U_D$, which influences the rescaling (7.1).

Another type of element in $\boldsymbol {A}$ are those that depend on velocity $v$ (and the background distribution function $H(v)$). They form continuous profiles. However, this continuity can be lost after the rescaling (7.1). Such elements enter the submatrices $\boldsymbol {C}^f$, $\tilde {\boldsymbol {F}}$ and $\boldsymbol {C}^E$, whose BE is described in §§ 7.5, 7.6 and 7.7, respectively.

7.2 The algorithm for BE $\hat {\boldsymbol {F}}$

7.2.1 Main idea

The main idea behind the algorithm for BE $\hat {\boldsymbol {F}}$ is to decompose this matrix into sets of elements such that all elements would have the same value within each set. After that, one extends large enough sets (as will be described further) in such a way that each extended set is encoded by a single STMC gate. However, after the extension, some sets may end up overlapping each other. In this case, the elements in the intersections (i.e. in the overlapping regions) should be corrected by additional STMC gates. The extension of some small sets that include a few matrix elements often leads to a significant overlapping with other sets. Therefore, it is often more efficient to encode a small set by encoding each element in it by its own STMC gate. Also, the representation of a matrix by sets is not unique. Because of that, we seek a minimal number of sets maximizing the number of large sets whilst minimizing the number of intersections and the number of small sets. Ideally, such optimization should be done using, for instance, deep reinforcement learning (see Sutton & Barto Reference Sutton and Barto2018). However, in the current version of the algorithm, the parameters are defined ad hoc.

This algorithm can be performed in $N_{\rm iter}$ iterations:

(7.2)\begin{equation} O_H \approx O_H^{(N_{\rm iter})}\ldots O_H^{(2)} O_H^{(1)}. \end{equation}

Here $O_H^{(1)}$ is the initial version of the oracle $O_H$, $O_H^{(k)}$ is the correction to $O^{(1)}_H$ provided by the $(k>1)$th iteration, and the final correction of the overlapping elements and encoding of small sets is performed at the $N_{\rm iter}$th iteration. After the $k$th iteration, one has

(7.3)\begin{equation} \delta \boldsymbol{A}^{(k)} = \boldsymbol{A} - \boldsymbol{A}^{(k)}, \end{equation}

where the matrix $\boldsymbol {A}^{(k)}$ is encoded by the sequence $O_H^{(k)}\ldots O_H^{(2)} O_H^{(1)}$. If $\delta A^{(k)}_{i_r, i_c} = 0$ then the first $k$ iterations correctly encode the matrix element at the row $i_r$ and the column $i_c$. If $|\delta A^{(k)}_{i_r, i_c}| > 0$ then the corresponding element should be corrected. To do that, one considers the matrix

(7.4)\begin{equation} \tilde{A}^{(k+1)}_{i_r, i_c} = \left\{\begin{array}{@{}ll} A^{(k)}_{i_r, i_c}, & |\delta A^{(k)}_{i_r, i_c}| > 0,\\ 0, & \delta A^{(k)}_{i_r, i_c} = 0, \end{array}\right. \end{equation}

which provides information about which matrix elements still should be corrected. At the $(k+1)$th iteration, one decomposes $\tilde {\boldsymbol {A}}^{(k+1)}$ into sets whose extension is restricted by the condition that the extended sets would not overlap the zero elements of $\tilde {\boldsymbol {A}}^{(k+1)}$. The $(k+1)$th iteration constructs the oracle $O_H^{(k+1)}$ correcting $\boldsymbol {A}^{(k)}$ by using $\tilde {\boldsymbol {A}}^{(k+1)}$. The correcting procedure is discussed in § 7.2.4. Because each extension of sets may lead to additional overlapping elements, it is not guaranteed that the $(k+1)$ iteration corrects all overlapping elements. For that reason, at the last iteration, the matrix $\tilde {\boldsymbol {A}}^{(N_{\rm iter})}$ is encoded element by element without extending the sets. In other words, the $N_{\rm iter}$th iteration corrects all remaining badly encoded matrix elements without introducing its own errors, which otherwise appear during the sets extension. Finally, by varying the number $N_{\rm iter}$, one can construct the oracle $O_H$ with a near-optimal scaling with respect to the matrix size.

7.2.2 Decomposition of the matrix

The elements in the submatrix $\hat {\boldsymbol {F}}$ are arranged along several diagonals (figure 1). The number of these diagonals is less or equal to the non-sparsity of the matrix $\boldsymbol {A}$ and depends neither on $N_x$ nor on $N_v$. Therefore, the first step is to represent $\hat {\boldsymbol {F}}$ as a group $\{\mathcal {D}\}$ of separate diagonals $\mathcal {D}$ (figure 12a) and then consider each diagonal $\mathcal {D}$ separately. One can specify a diagonal $\mathcal {D}$ by using a constant integer $\Delta i$ as

(7.5)\begin{equation} \mathcal{D} = \{\hat{F}_{i_r, i_r + \Delta i}\,|\, i_r = [\delta_{{\rm sign}_{\Delta i}, -1}|\Delta i|, N_{xv} - \delta_{{\rm sign}_{\Delta i}, 1}\Delta i)\}, \end{equation}

where ${\rm sign}_{\Delta i} = \Delta i / |\Delta i|$. Each $\mathcal {D}$ is encoded independently of other diagonals (as shown by various dashed blocks in figure 13) that results in a linear dependence of the depth of the oracle $O_H$ on the non-sparsity $\varsigma _{\hat {\boldsymbol {F}}}$ of the submatrix $\hat {\boldsymbol {F}}$.

Figure 12. A schematic showing decomposition of a matrix into sets for the construction of the oracle $O_H$. (a) A simple matrix of size $N_A = 8$ is taken as an example. This matrix consists of three diagonals marked in different colours, and each diagonal is considered separately. The element values are indicated as $e_i$, where $e_i \neq e_j$ if $i\neq j$. (b) A schematic showing how the main diagonal marked in red in (a) can be split in several sets where each set groups several matrix elements of the same value. (c-i) A matrix with bitstrings of the row indices of the original set $S_2$ and of its extended version $S^{\rm ext}_{2}$. This extension results in three overlapping elements marked in orange, i.e. $S^{\rm ext}_2$ overlaps elements of the sets $S_{0,0}$, $S_1$ and $S_{0,1}$. (c-ii) A bitstring matrix of the set $S_2$ that has been split into two sets, $S_{2,0}$ and $S_{2,1}$, where only the latter is extended. This extension results in a single overlapping element. (d) A schematic showing the splitting and extension (if necessary) of sets in the left and right diagonals. The empty cells indicate that the considered diagonal does not have elements at the corresponding rows. The bits indicated in grey in the bitstring matrices are chosen as the control nodes of the STMC gates encoding the extended sets (figure 13).

Figure 13. The circuit of the oracle $O_H$ encoding the elements of the tridiagonal matrix from figure 12(a). The dashed blocks indicate the parts of the oracle encoding the main diagonal (red box), the left diagonal (green box) and the right diagonal (blue box), correspondingly. Here, the register $r$ encodes the matrix row indices. The ancilla register $a_r$ has a similar meaning to that explained in (6.25) (although here the register has two qubits) and is used to address the matrix diagonals. The ancilla $a_e$ is initialized in the zero state, and the matrix elements’ values are written into the amplitude of $\mid {0}\rangle_{a_e}$. The gate $[e_i]$ is a schematic representation of a rotation gate whose rotation angle is chosen to encode the value $e_i$ as explained in (6.13). The gate $[e_i^*]$ is a rotation gate whose angle is computed taking into account the overlapping of the element $e_i$ with one or several extended sets. For instance, the third gate in the red box encodes the set $S_{2,0}$. The fourth gate encodes the extended set $S_{2,1}^{\rm ext}$. Since $S^{\rm ext}_{2,1}$ intersects with $S_{0,1}$, the angle for the fifth gate is computed taking into account the action of the fourth gate, as explained in § 7.2.4.

Each diagonal is represented by a group of sets $\{S\}$ where each set $S$ contains elements of the same value $v^{(S)}$, i.e.

(7.6)\begin{equation} S = \{v^{(S)}, \{i_r\}\,|\, \hat{\boldsymbol{F}}_{k, k + \Delta i} = v^{(S)},\ \forall k \in \{i_r\} \}, \end{equation}

where $\{i_r\}$ are the row indices of all matrix elements described by the set $S$. Each set stores only $v^{(S)}$ and $\{i_r\}$. A set can be divided into smaller sets (§ 7.2.5) and, thus, two different $S$ can have the same $v^{(S)}$.

Once the matrix $\hat {\boldsymbol {F}}$ is split into $\{\mathcal {D}\}$, where $\mathcal {D} = \{S\}$, one constructs a circuit representation of $O_H$. The element-by-element encoding is inefficient, because it will require $\mathcal {O}(N_{xv})$ quantum gates. The main purpose of our algorithm is to construct a circuit with the scaling better than $\mathcal {O}(N_{xv})$. Ideally, the scaling should be $\mathcal {O}(\text {polylog}(N_{xv})\varsigma _{\hat {\boldsymbol {F}}})$. (As shown in Novikau et al. (Reference Novikau, Dodin and Startsev2023), the QSVT for solving a multi-dimensional stationary wave problem described by a matrix of size $N$ and non-sparsity $\varsigma$ can have a polynomial speedup in comparison to classical conjugate-gradient-based algorithms if the corresponding BE oracle scales as $\mathcal {O}(\text {polylog}(N)\varsigma )$).

To achieve a better scaling, the algorithm encodes the sets instead of the matrix elements by using the fact that each $S$ contains elements with the same value. The simplest way to encode $S$ is to use a quantum arithmetic comparator (see Suau, Staffelbach & Calandra Reference Suau, Staffelbach and Calandra2021; Novikau et al. Reference Novikau, Dodin and Startsev2023) to find whether $i_r \in \{i_r\}_{S}$, where $i_r$ is the row index encoded into the registers $r_v$, $r_x$ and $r_f$ (§ 4.1). However, for each $S$, one will need to use at least two comparators, and each comparator requires two ancilla qubits (although the ancillae can be potentially reused by various comparators). Instead, we try to arrange sets in such a manner that each large enough set would be encoded by a single STMC gate. Such an arrangement is performed by extending $S$. (Another possibility could be to use quantum arithmetic operators to encode the binary representation of the matrix elements. However, that would require much more qubits and be hard to test numerically, which is why we do not use this approach here.)

7.2.3 Set extension

To encode all elements of $S$ at once, the algorithm extends $S$ in such a way that the extended set $S^{\rm ext}$ is computed by a single STMC gate. To make it possible, it is necessary to complete $\{i_r\}$ of $S$ by another set $\{i_r^{\rm com}\}$ such that the resulting set $\{i_r\}\cup \{i_r^{\rm com}\}$ could be represented by (6.20) and (6.21). In other words, $S$ is extended to ensure that it would be possible to find a STMC gate $C_{\{q_{c\delta }\}}G^{(a_e)}$ such that $G_{00} = v^{(S)}$ and $\mathcal {S} \equiv \{i_r\}\cup \{i_r^{\rm com}\}$.

The extended set $S^{\rm ext}$ consists of a core, which is the original non-extended set $S$, and a complement $S^{\rm com}$: $S^{\rm ext} = S \cup S^{\rm com}$ and $S \cap S^{\rm com} = \varnothing$, where $S^{\rm com}$ includes $N_c$ elements described by the row indices $\{i_r^{\rm com}\}$:

(7.7)\begin{equation} S^{\rm ext} = \{v^{(S)}, \{i_r\}\cup\{i_r^{\rm com}\}\,|\, \hat{F}_{k, k + \Delta i} = v^{(S)},\ \forall k \in \{i_r\} \}. \end{equation}

Note here that the matrix elements $\hat {F}_{j, j + \Delta i}$ for $j \in \{i_r^{\rm com}\}$ may have a value different from $v^{(S)}$ but the extended set $S^{\rm ext}$ assumes that they do have the value $v^{(S)}$. Thus, the extended sets can assign wrong values to some matrix elements. In this case, we speak about the overlapping of several sets. The cores of extended sets never intersect, $S_j \cap S_k = \varnothing$ for $j\neq k$, but the complements can overlap with the cores or complements of other sets. Hence, one can have only complement–complement and complement–core overlapping. The matrix elements that sit in the overlapping regions of several sets should be corrected. The correction employs supplemental STMC gates. Therefore, it is important to minimize the number of overlapping elements. The correction can be done in different ways and a possible algorithm for that is described in § 7.2.4.

We should also note that, sometimes, small sets are extended to significantly larger sets, thus, drastically increasing the number of overlapping elements. Because of that, it is better to encode small enough sets such that each element in them is computed by a separate STMC gate.

In well-structured matrices, such as those that appear in classical linear wave problems, the number of sets is significantly less than the number of matrix elements. Since each extended set is encoded by a single STMC gate, the scaling of the resulting oracle should be significantly better than in the case when the matrix is encoded element by element. For instance, let us assume that a diagonal $\mathcal {D}$ consists of $N$ elements. If all these elements are equal to each other then $\mathcal {D}$ comprises a single set, and this set can be encoded into a quantum circuit by using a single STMC rotation gate. Another case is when there are $N_S$ elements of the same value where $(N - N_S) \ll N$, then $\mathcal {D}$ has a dominant set with $N_S$ elements. This set can be extended to the whole diagonal as shown in figure 12(b) for the set $S_2$. In this case, the complement $S^{\rm com}$ contains $(N-N_S) \ll N$ elements and each of these elements are encoded by a separate STMC gate. This means that the oracle $O_H$ encoding the diagonal $\mathcal {D}$ has $(N-N_S + 1) \ll N$ STMC gates.

However, extending the original set may not always be the best approach. Instead, $S$ can be split (§ 7.2.5) into several subsets and then each subset can be extended individually. As demonstrated in figure 12(b) for the set $S_2$, the splitting allows us to reduce the number of overlapping elements from $N_c = 3$ to $N_c = 1$.

7.2.4 Correction of overlapping elements

If two sets intersect with each other, the overlapping elements should be corrected by supplemental gates. Let us consider the correction on the example from figure 12, where, for instance, two sets, $S_{2,1}^{\rm ext}$ and $S_{0,1}$, overlap at the element $e_0$. According to figure 12(c-ii), to encode the set $S_{2,1}^{\rm ext}$, one needs a single STMC gate controlled by the most significant qubit of the input register $r$ (as shown in figure 13 by the fourth gate in the red box). To encode an arbitrary complex value, one can use the gate $R_c(2\theta _{z,2}, 2\theta _{y,2})$ described in Appendix C, which should act at the zero state of the ancilla $a_e$. In particular, to compute the value $e_2$, one needs to chose the gates’ parameters in such a way that $e_2 = \cos (\theta _{y,2})\exp ({-\mathrm {i}\theta _{z,2}})$. In this case, one obtains

(7.8)\begin{equation} \mid{\psi_{\rm ext}}\rangle_{a_e} = R_c(2\theta_{z,2}, 2\theta_{y,2})\mid{0}\rangle_{a_e} = e_2 \mid{0}\rangle_{a_e} + w_2 \mid{1}\rangle_{a_e}, \end{equation}

where $w_2$ is a complex value defined by the angles of the gate $R_c$. Thus, two components of the state $\psi _{\rm ext}$ are $\psi _{\rm ext,0} = e_2$ and $\psi _{\rm ext,1} = w_2$. This gate controlled by the most significant qubit of the register $r$ (figure 13) entangles the state $\mid {\psi _{\rm ext}}\rangle_{a_e}$ with the states $\mid {i_r}\rangle_{r}$, where $i_r = [4,7]$. However, at $i_r = 7$, the main diagonal has an element with the value $e_0$ instead of $e_2$. This means that a supplemental gate $R_c(2\theta _{z,c}, 2\theta _{y,c})$ should act on the state $\mid {\psi _{\rm ext}}\rangle_{a_e}\mid {7}\rangle_{r}$ to encode the value $e_0$:

(7.9)\begin{equation} \mid{\psi_{\rm res}}\rangle_{a_e} = R_c(2\theta_{z,c}, 2\theta_{y,c}) \mid{\psi_{\rm ext}}\rangle_{a_e} = e_0 \mid{0}\rangle_{a_e} + w_0 \mid{1}\rangle_{a_e}. \end{equation}

Here we are not interested in the value $w_0$. This explains the last gate in the red box in figure 13, denoted $[e_0^*]$.

Equation (7.9) leads to the following equations for $\theta _{z,c}$ and $\theta _{y,c}$:

(7.10a)\begin{gather} g_{00}\cos\theta_{y,c} - g_{01}\sin\theta_{y,c} = \psi_{\rm res,0}^{\rm real}, \end{gather}

(7.10b)\begin{gather}g_{10}\cos\theta_{y,c} - g_{11}\sin\theta_{y,c} = \psi_{\rm res,0}^{\rm imag}. \end{gather}

Here

(7.11a)\begin{gather} g_{00} = \psi_{\rm ext,0}^{\rm real}\cos\theta_{z,c} + \psi_{\rm ext,0}^{\rm imag}\sin\theta_{z,c}, \end{gather}

(7.11b)\begin{gather}g_{01} = \psi_{\rm ext,1}^{\rm real}\cos\theta_{z,c} - \psi_{\rm ext,1}^{\rm imag}\sin\theta_{z,c}, \end{gather}

(7.11c)\begin{gather}g_{10} = \psi_{\rm ext,0}^{\rm imag}\cos\theta_{z,c} - \psi_{\rm ext,0}^{\rm real}\sin\theta_{z,c}, \end{gather}

(7.11d)\begin{gather}g_{11} = \psi_{\rm ext,1}^{\rm real}\sin\theta_{z,c} + \psi_{\rm ext,1}^{\rm imag}\cos\theta_{z,c}. \end{gather}

In the more general case, when the $N_g$ gates act on the same element $e_{\rm res}$, then one can find $\theta _{y,c}$ and $\theta _{z,c}$ in (7.10)–(7.11) from

(7.12a)\begin{gather} \psi_{\rm ext} = R_{c,N_g-1} R_{c, N_g-2}\cdots R_{c,1} R_{c,0}\mid{0}\rangle_{a_e}, \end{gather}

(7.12b)\begin{gather}\psi_{\rm res, 0} = e_{\rm res}. \end{gather}

7.3 Extrapolation of the matrix ${D}^A$

To construct the oracle $O_H$, one needs to compute the matrix $\boldsymbol {D}^A$ related to the oracle $U_D$ (6.9) to perform the rescaling (6.13). The parts of $\boldsymbol {D}^A$ related to the submatrices $\boldsymbol {C}^E$, $\boldsymbol {C}^f$ and $\boldsymbol {S}$ have a trivial structure and can be predicted for any $N_x$ and $N_v$. For instance, $D^A_{N_{xv} + i_x N_v + i_v,\ N_{xv}+ i_x N_v + i_v} = 1$ for $i_x = [0, N_x)$ and $i_v = [1, N_v)$. These elements correspond to the main diagonal of the submatrix $\boldsymbol {S}$, where the matrix non-sparsity is $1$. In the rows with indices $N_{xv} + i_x N_v$, the column indices of the non-zero elements are computed by a single Hadamard gate acting on the ancilla $a_f$ and by $n_v$ Hadamard gates acting on all qubits of the register $a_v$. (These Hadamard gates appear in both $O_F$ and $O_F^{\dagger}$, according to (6.9) and figure 10.) Because of that, $D^A_{N_{xv} + i_x N_v,\ N_{xv}+ i_x N_v} = 1/2$ that corresponds to the diagonal elements at $i_v = 0$ in the submatrix $\boldsymbol {S}$, and $D^A_{N_{xv} + i_x N_v,\ i_x N_v + i_v} = 1/2^{n_v/2+1}$ with $i_v = [0, N_v)$ that corresponds to the non-zero elements of the submatrix $\boldsymbol {C}^f$.

Let us denote the upper left $N_{xv}\times N_{xv}$ part of $\boldsymbol {D}^A$ as $\boldsymbol {D}^F$. The submatrix $\boldsymbol {D}^F$ contains elements $d_{i_r i_c}$ by which the submatrix $\boldsymbol {F}$ should be rescaled in (6.13). The computation of $\boldsymbol {D}^F$ is numerically challenging for a large $N_{xv}$. Instead, $\boldsymbol {D}^F$ of a large size can be extrapolated by using known matrices $\boldsymbol {D}^F$ of several small sizes. To do that, one creates a template of $\boldsymbol {D}^F$ by using computed $\boldsymbol {D}^F$ for small $N_x$ and $N_v$. This can be done because the relative positions of the submatrix elements and the elements’ values do not change with the increase of $N_{xv}$. The template is used to reconstruct $\boldsymbol {D}^F$ for any $N_{xv}$. The algorithm for the creation of this template is the following.

Let us consider the matrix $\boldsymbol {D}^F$ that is split into several diagonals in the same way as described in § 7.2.2. The algorithm considers each of these diagonals independently and assumes that the row index $r$ of this matrix depends on $N_{\rm ind}$ indices $r_l$:

(7.13)\begin{equation} r = \sum_{l=0}^{N_{\rm ind}-1} r_l \prod_{j=(l+1)}^{N_{\rm ind}} N_j. \end{equation}

Here each index $r_l = [0, N_l)$ corresponds to the $l$th dimension of the size $N_l$ and $N_{N_{\rm ind}} = 1$. For instance, any matrix describing dynamics in three-dimensional real space ($N_{\rm ind} = 3$) has $r = r_0 N_y N_x + r_1 N_x + r_2$. In our case, $r_0 = i_{xr}$ and $r_1 = i_{vr}$ according to (6.24). Also, we define the following ordered set of indices:

(7.14)\begin{equation} \overrightarrow{\{r\}}_{l} = \{r_j \,|\, r_j = [0, N_j),\ j = [0, l) \}. \end{equation}

The ordered set $\overrightarrow {\{r\}}_{l}$ is empty if $l = 0$.

Let us consider the diagonal $\mathcal {D}$ defined in (7.5), where instead of $\hat {\boldsymbol {F}}$ we now take $\boldsymbol {D}^F$. The vector with all matrix elements from $\mathcal {D}$ is denoted $\boldsymbol {M}$. It has $N_{xv}$ elements, where we formally define $M_{r} = \varnothing$ if the diagonal does not have any element at the row $r$. Due to the dependence (7.13), $\boldsymbol {M}$ has a nested structure with $N_0$ blocks, where each block has $\prod _{j=1}^{N_{\rm ind}} N_j$ rows. These rows can be combined into $N_1$ sub-blocks, where each sub-block has $\prod _{j=2}^{N_{\rm ind}} N_j$ rows, and so on. Let us call all blocks of an equal size a layer. There are $N_{\rm ind}$ layers. The $i_l$th block in the $l$th layer is defined as

(7.15)\begin{equation} B_{\overrightarrow{\{r\}}_{l}, i_l} = \{B_{\overrightarrow{\{r\}}_{l} \cup i_l, i_{l+1}}\}, \end{equation}

where $i_l = [0, N_l)$, $l = [0, N_{\rm ind})$ and each block of the $(N_{\rm ind} - 1)$th layer is a single matrix element, i.e.

(7.16)\begin{equation} B_{\overrightarrow{\{r\}}_{q}, i_q} = \left\{M_{r} \,|\, r = \sum_{r_l\in\overrightarrow{\{r\}}_{q}} r_l \prod_{j=(l+1)}^{q} N_j + i_q\right\}, \end{equation}

where $q = N_{\rm ind} - 1$. Hence, the $l$th layer has $N_l$ blocks, and each of these blocks has $\prod _{j=l+1}^{N_{\rm ind}} N_j$ matrix elements.

Identical adjacent blocks in the $l$th layer are combined into blocksets. A blockset in the $l$th layer is denoted $\mathcal {B}_l$ and is defined as

(7.17)\begin{equation} \mathcal{B}_l = \{B_{\overrightarrow{\{r\}}_{l}, r_b}, r_b, r_e\,|\, B_{\overrightarrow{\{r\}}_{l}, r} = B_{\overrightarrow{\{r\}}_{l}, r_b},\ \forall r = [r_b, r_e)\}, \end{equation}

where each blockset in the $q = (N_{\rm ind} - 1)$th layer combines identical neighbour matrix elements

(7.18)\begin{equation} \mathcal{B}_q = \{M_{r}, r_b, r_e\,|\, M_{r} = M_{r_b}, \forall r = [r_b, r_e)\}. \end{equation}

Since each block in the $l$th layer includes blocks from the $(l+1)$th layer (7.15), each block can be described by a group of blocksets, i.e.

(7.19)\begin{equation} B_{\overrightarrow{\{r\}}_{l},i_l} = \{\mathcal{B}_{l+1}\}_{i_l}, \end{equation}

where we keep the index $i_l$ to indicate that $\{\mathcal {B}_{l+1}\}_{i_l}$ represents the $i_l$th block at the $l$th layer. If the block $B_{\overrightarrow {\{r\}}_{-1},i_{-1}}$ is formally defined as the whole diagonal $\boldsymbol {M}$, then $\{\mathcal {B}_{0}\}_{i_{-1}}$ is a group of blocksets where the blocks of the zeroth layer are sorted. By combining (7.17) and (7.19), one can see that each blockset in the $l$th layer can be defined as a group of blocksets from the $(l+1)$th layer:

(7.20)\begin{equation} \mathcal{B}_l = \{\{\mathcal{B}_{l+1}\}_{r_b}, r_b, r_e\,|\, \{\mathcal{B}_{l+1}\}_{r} = \{\mathcal{B}_{l+1}\}_{r_b},\ \forall r = [r_b, r_e)\}. \end{equation}

Thus, the whole diagonal $\boldsymbol {M}$ can be represented as nested blocksets. Representing a matrix by a group of diagonals comprising blocksets can be regarded as creating a template that represents a compressed image of the matrix.

In the matrix $\boldsymbol {D}^F$, the number of diagonals $\mathcal {D}$ does not change with $N_{xv}$. The matrix elements $M_r$ stored in the blocksets (7.18) are also independent of $N_{xv}$. (Some matrix elements of $\boldsymbol {D}^A$ do change with $N_v$ because the non-sparsity of $\boldsymbol {C}^f$ and $\boldsymbol {C}^E$ changes with $N_v$. Yet, the non-sparsity of the submatrix $\boldsymbol {F}$ does not change with $N_{xv}$, and because of that, the elements of $\boldsymbol {D}^F$ do not change as well.) The indices $r_b$ and $r_e$ in (7.20) and (7.18) depend linearly on $N_l$, i.e.

(7.21a)\begin{gather} r_b = \alpha_b + \beta_b N_l, \end{gather}

(7.21b)\begin{gather}r_e = \alpha_e + \beta_e N_l, \end{gather}

where $l=[0,N_{\rm ind})$. The unknown coefficients $\alpha _b$, $\alpha _e$, $\beta _b$ and $\beta _e$ are different for different blocksets. To compute these coefficients, one constructs the templates for $N_{\rm ind}+1$ matrices where each matrix has at least one $N_l$ different from the corresponding $N_l$ of all other matrices. (If needed, each of the $N_{\rm ind}+1$ matrices can have all $N_{\rm ind}$ dimensions different from the dimensions of other matrices.) For each of these computed $(N_{\rm ind}+1)$ templates, one knows the row indices $r_b$ and $r_e$ for each blockset. One substitutes these indices to (7.21) from which the unknown coefficients can now be calculated.

Once the coefficients are computed, one can construct $\boldsymbol {D}^F$ of an arbitrary size $N_{xv}$. Equations (7.21) allow us to compute the positions $r_b$ and $r_e$ for all blocksets in $\boldsymbol {D}^F$ for the required $N_{xv}$.

For the matrix $\boldsymbol {D}^F$, $N_{\rm ind}$ equals two, where $N_0 = N_x$ an $N_1 = N_v$. Thus, one needs to precalculate three matrices $\boldsymbol {D}^F$ with various $N_x$ and $N_v$ to use (7.21). Yet, one should keep in mind that $N_x$ and $N_v$ should be large enough to take into account the influence of the boundary elements in $\boldsymbol {D}^F$ properly. More precisely, one should consider only matrices with $n_x \geq 4$ and $n_v \geq 4$.

As an example illustrating the decomposition of a matrix into blocksets, let us consider a diagonal matrix $\boldsymbol {D}^{\rm ex}$ with diagonal elements $D^{\rm ex}_{i, i}$, where $i = i_y N_x + i_x$ with $i_y = [0,N_y)$ and $i_x = [0,N_x)$ for $N_y = 2^{n_y}$ and $N_x = 2^{n_x}$:

(7.22a)\begin{gather} D^{\rm ex}_{i, i} = e_0,\quad i_y = [0,N_y/2),\quad i_x = [0,N_x-2], \end{gather}

(7.22b)\begin{gather}D^{\rm ex}_{i, i} = e_1,\quad i_y = [0,N_y/2),\quad i_x = N_x-1, \end{gather}

(7.22c)\begin{gather}D^{\rm ex}_{i, i} = e_0,\quad i_y = [N_y/2,N_y),\quad i_x = 0\text{ and } N_x-1, \end{gather}

(7.22d)\begin{gather}D^{\rm ex}_{i, i} = e_1,\quad i_y = [N_y/2,N_y),\quad i_x = [1,N_x-2]. \end{gather}

Here, $e_0$ and $e_1$ are some non-equal values. Due to the dependence on the indices $i_x$ and $i_y$, the main diagonal has a nested structure with $N_y$ blocks and with $N_x$ elements in each block. Other diagonals are empty (the matrix elements in these diagonals are equal to zero). The main diagonal for different $n_y$ and $n_x$ is shown in figure 14.

Figure 14. A schematic showing elements in the main diagonal of the matrix $\boldsymbol {D}^{\rm ex}$ described in (7.22) for various $N_x = 2^{n_x}$ and $N_y = 2^{n_y}$. The matrix has $N_y$ blocks with $N_x$ elements each. These blocks are grouped into two blocksets, $\mathcal {B}_0^I$ and $\mathcal {B}_0^{II}$. The elements in each block of the blockset $\mathcal {B}_0^I$ are combined into two blocksets indicated by different shades of green. The elements in each block of the blockset $\mathcal {B}_0^{II}$ are combined into three blocksets indicated by different shades of blue.

In the zeroth layer the diagonal is split into two blocksets, $\mathcal {B}_0^I$ and $\mathcal {B}_0^{II}$, independently of the matrix size. In the next layer, $\mathcal {B}_0^I$ contains two blocksets (which are indicated by different shades of green in figure 14), and $\mathcal {B}_0^{II}$ contains three blocksets (indicated by different shadows of blue). Although the number of blocks and elements in $\boldsymbol {D}^{\rm ex}$ changes with $N_x$ and $N_y$, the number of blocksets in any layer does not. For instance, independently of $N_y$, the blockset $\mathcal {B}_0^{I}$ always contains only two blocksets. Indeed, according to (7.17), $\mathcal {B}_0^{I}$ stores a single copy of one of $N_y/2$ identical blocks in the first half of the diagonal. In turn, this single copy contains two blocksets, one of which saves a single copy of the elements $e_0$ and another saves the element $e_1$.

The indices stored by the blockset $\mathcal {B}_0^{I}$ are

(7.23a,b)\begin{equation} r_b = 0,\quad r_e = \tfrac{1}{2} N_y. \end{equation}

The indices of the first inner blockset of $\mathcal {B}_0^{I}$ are $r_b = 0$ and $r_e = N_x-2$. The indices of the second inner blockset are $r_b = N_x-2$ and $r_e = N_x-1$.

The blockset $\mathcal {B}_0^{II}$ have the following boundaries:

(7.24a,b)\begin{equation} r_b = \tfrac{1}{2} N_y,\quad r_e = N_y. \end{equation}

The indices of the inner blocksets of $\mathcal {B}_0^{II}$ can be easily found from figure 14 or (7.22).

7.4 Encoding the submatrix ${S}$

To encode the submatrix $\boldsymbol {S}$ described in (4.12), we use the two gates $R_x$ with the following angles:

(7.25a)\begin{gather} \theta_{S,b} ={-}2 \arcsin(\omega_0), \end{gather}

(7.25b)\begin{gather}\theta_{S,e} ={-}2 \arcsin(2 \omega_0) - \theta_{S,b}. \end{gather}

As one can see from the red dashed box in figure 9, the elements of $\boldsymbol {S}$ are encoded by three gates, $R_x(\theta _{S,b})$, $R_x(\theta _{S,e})$ and $X$, applied to the ancilla $a_e$. The first gate $R_x(\theta _{S,b})$ encodes $\mathrm {i}\omega _0$ into the amplitude of the zero state of the ancilla $a_e$. Due to the rescaling (7.1), the value $\mathrm {i}\omega _0$ of each element at $i_{vr} = 0$ in $\boldsymbol {S}$ is multiplied by the factor $2$. This is taken into account by the second gate, $R_x(\theta _{S,e})$, which corrects the action of $R_x(\theta _{S,b})$.

7.5 Encoding the submatrix ${C}^f$

To encode the matrix elements of $\boldsymbol {C}^f$ (4.11) that depend on the velocity, we use the circuit shown in figure 15. It encodes the function $\sin (\phi _i)$ for $i = [0, N_t)$ with $N_t = 2^{n_t}$:

(7.26)\begin{equation} \phi_i = \alpha_0 + {\rm i} \Delta\phi. \end{equation}

Here $n_t$ is the number of qubits in the register $r_t$ and $\Delta \phi = 2\alpha _1/N_t$. This circuit scales as $\mathcal {O}(n_t)$.

Figure 15. The circuit encoding $\sin (\phi _i)$ according to (7.26). Here, $R_k = R_y(2\alpha _1/2^k)$.

To encode the velocity grid $v_i$, where $v_{\rm max}\ll 1$ due to the normalization (6.7), one can use the circuit 15, since $\sin (\phi _i) = \phi _i + \mathcal {O}(\phi _i^3)$. In this case, $n_t = n_v$, $\alpha _0 = - v_{\rm max}$ and $\alpha _1 = |\alpha _0| N_v/ (N_v - 1)$. All non-zero elements of the submatrix $\boldsymbol {C}^f$ are encoded by using a single call to the circuit 15.

7.6 Encoding the submatrix $\tilde {{F}}$

The elements of $\tilde {\boldsymbol {F}}$ described after (4.9) depend on $v$ linearly and are encoded using the same technique described in § 7.5.

7.7 Encoding the submatrix ${C}^E$

To encode the matrix elements of $\boldsymbol {C}^E$ (4.10), one can use the QSVT to approximate the odd function $v H(v)$. First of all, the odd function $v H(v)$ is approximated by the form

(7.27)\begin{equation} f_{vH}(y_i) = p_1 \arcsin(y_i) \exp({p_0\arcsin(y_i)^2}), \end{equation}

where $y_i = \sin (\phi _i)$ for $i = [0, N_v)$ where the angles $\phi _i$ are computed using (7.26) with the parameters $\alpha _0 = -1$ and $\alpha _1 = |\alpha _0| N_v/(N_v-1)$. The parameter $p_0$ depends neither on $n_x$ nor on $n_v$. The parameter $p_1$ decreases linearly with $N_v$ and does not depend on $N_x$. The approximation $f_{vH}(y_i) = v_i H(v_i) \pm \varepsilon _{\rm ce, appr}$ is found by solving a nonlinear least-square problem.

The function $f_{vH}(y_i)$ is computed by using QSVT, where the variable $y_i$ is encoded by the circuit 15. Because flat background temperature and density profiles are considered, this function does not depend on $x$ and can be encoded for all $x$ points by a single call to the QSVT circuit. The absolute error $\epsilon _{\rm CE, qsvt}$ of the QSVT computation of the function (7.27) rapidly decreases with the number of QSVT angles $N_{\rm CE,qsvt}$, e.g. $\epsilon _{\rm CE, qsvt} = 10^{-6}$ for $N_{\rm CE,qsvt} = 16$. Each QSVT angle is associated with a single call to the oracle 15 that scales linearly with $n_v$. Thus, the scaling of the circuit encoding the submatrix $\boldsymbol {C}^E$ is $\mathcal {O}(n_v \log _2(\epsilon _{\rm CE, qsvt}^{-1}))$.

If the temperature and density depend on the spatial coordinate, then one needs to implement additional QSVT circuits that would encode the change in the amplitude and the width of the bumps of the function $v H$. Another option is to substitute these QSVT circuits by a single one by applying multivariable QSP (see Rossi & Chuang Reference Rossi and Chuang2022) that operates with several BE oracles at once, thus computing a polynomial of multiple variables, i.e. $x$ and $v$ in our case.

7.8 Scaling of the BE oracle

The final circuit of the BE oracle encoding the matrix $\boldsymbol {A}$ can be found in Novikau (Reference Novikau2024b) and was numerically tested using the QuCF framework (Novikau Reference Novikau2024c). Now, let us estimate the scaling of this circuit. It comprises the scaling of the oracle $U_D$ and the oracle $O_H$. As discussed in § 6.6, $U_D$ scales as $\mathcal {O}(n_x + n_v)$ due to the arithmetic operators. The circuit of $O_H$ consists of several pieces: the oracle encoding the submatrix $\hat {\boldsymbol {F}}$, whose encoding is performed by the procedure discussed in § 7.2, and the oracles encoding the submatrices $\boldsymbol {S}$ (§ 7.4), $\boldsymbol {C}^f$ (§ 7.5), $\tilde {\boldsymbol {F}}$ (§ 7.6) and $\boldsymbol {C}^E$ (§ 7.7). The scaling of the circuit for $\hat {\boldsymbol {F}}$ is shown in figure 16 and is estimated as $\mathcal {O}(\varsigma _{\hat {\boldsymbol {F}}}{\rm poly}(N_x) {\rm poly}(n_v))$. The poor scaling with respect to $N_x$ is caused by the matrix elements that appear due to the non-zero diffusivity $\eta$ (2.11). The non-sparsity $\varsigma _{\hat {\boldsymbol {F}}}$ of the submatrix $\hat {\boldsymbol {F}}$ does not change with $N_x$ or $N_v$. Yet, $\varsigma _{\hat {\boldsymbol {F}}}$ depends on the discretization method or boundary conditions of the simulated kinetic model.

Figure 16. The dependence of the number of non-zero elements $N_{\rm nz}$ in the matrix $\boldsymbol {F}_F$ on (a) $n_v$ for various $n_x$ and on (b) $n_x$ for various $n_v$. (c) The dependence of the number of STMC gates in the oracle $O_H$ necessary for encoding $\boldsymbol {F}_F$ on $n_v$ for various $n_x$. (d) The dependence of the number of STMC gates on $n_x$ for various $n_v$. Here, the text in different colours indicate the fitting equations approximating the scaling with respect to $n_v$ and $n_x$ or $N_v$ and $N_x$.

The depth of the circuit for $\boldsymbol {S}$ scales as $\mathcal {O}(1)$ if one uses STMC gates, and, according to Claudon et al. (Reference Claudon, Zylberman, Feniou, Debbasch, Peruzzo and Piquemal2023), an arbitrary STMC gate controlled by $n$ qubits can be decomposed into a circuit with $\mathcal {O}(\log _2(n)^{\log _2(12)}\log _2(1/\epsilon _{\rm STMC}))$ depth, as mentioned in § 6.4. The oracles for $\boldsymbol {C}^f$ and $\tilde {\boldsymbol {F}}$ scale as $\mathcal {O}(n_v)$. The scaling of the oracle $\boldsymbol {C}^E$ implemented using QSVT is $\mathcal {O}(n_v\log _2(\epsilon _{\rm CE, qsvt}^{-1}))$.

Then, in summary, the main contribution to the scaling of the BE oracle of the whole matrix $\boldsymbol {A}$ is due to the oracle encoding $\hat {\boldsymbol {F}}$, and the resulting scaling is

(7.28)\begin{equation} \mathcal{O}[n_x + n_v + \varsigma_{\hat{\boldsymbol{F}}}{\rm poly}(N_x) {\rm poly}(n_v) + n_v\log_2(\epsilon_{\rm CE, qsvt}^{{-}1})]. \end{equation}