2.1. Semidiscretized central difference scheme

Let us consider a six-dimensional phase space $\Omega := [0, L] \times [0, L] \times [0, L] \times [-V, V] \times [-V, V] \times [-V, V]$ . Here, $L, V \in \mathbb{R}_+$ , and for simplicity, we assume that the three-dimensional (3-D) physical space is uniform in the range $[0, L]$ and the 3-D velocity space is uniform in the range $[-V, V]$ . Each direction of the domain is discretized with $N$ uniformly distributed points. Thus, the grid intervals are $\Delta x = \Delta y = \Delta z = L / N$ and $\Delta v := \Delta v_x = \Delta v_y = \Delta v_z = 2V / N$ where $N = 2^n$ and $n \in \mathbb{N}$ .

We define the six-dimensional scalar $f({\boldsymbol{x}},{\boldsymbol{v}})$ in the phase space $\Omega$ . We hereafter label the grid points in $\Omega$ with vectors of indices,

(2.1) \begin{equation} {\boldsymbol{j}}:=(j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z})\in \mathcal{I}_6:=[N]^{\otimes 6}, \end{equation}

where, $j_x, j_y$ and $j_z$ represent the index of the grid points in the 3-D physical space, and $j_{v_x}, j_{v_y}$ and $j_{v_z}$ represent the indices of the grid points in the 3-D velocity space, with each index taking a value from $0$ to $N-1$ . For convenience in later use, we define a flattened index for the multidimensional information as

(2.2) \begin{equation} j:=j_x+j_yN+j_zN^2+j_{v_x}N^3+j_{v_y}N^4+j_{v_z}N^5. \end{equation}

Now, the Vlasov equation for collisionless plasma is represented in a conservative form with the plasma distribution function in phase space as $f({\boldsymbol{x}},{\boldsymbol{v}},t)$ ,

(2.3) \begin{equation}\frac{\partial f({\boldsymbol{x}}, {\boldsymbol{v}}, t)}{\partial t}+ {\boldsymbol{v}} \boldsymbol{\cdot} \frac{\partial f({\boldsymbol{x}}, {\boldsymbol{v}}, t)}{\partial {\boldsymbol{x}}}+ \frac{q}{m} \left( {\boldsymbol{E}}({\boldsymbol{x}}, t) + {\boldsymbol{v}} \times {\boldsymbol{B}}({\boldsymbol{x}}, t) \right) \boldsymbol{\cdot} \frac{\partial f({\boldsymbol{x}}, {\boldsymbol{v}}, t)}{\partial {\boldsymbol{v}}}= 0.\end{equation}

Here, ${\boldsymbol{E}}$ , ${\boldsymbol{B}}$ and ${\boldsymbol{v}}$ are 3-D vectors, representing the electric field, magnetic field and velocity in phase space, respectively. To discretize the spatial derivative terms, we use the central difference method. For example, the spatial derivative along the $x$ direction can be written with second-order accuracy as

(2.4) \begin{equation}\frac{\partial f({\boldsymbol{x}}, {\boldsymbol{v}})}{\partial {\boldsymbol{x}}}= \frac{1}{2\Delta x} \left( f({\boldsymbol{x}} + \Delta x {\boldsymbol{e}}, {\boldsymbol{v}}) - f({\boldsymbol{x}} - \Delta x {\boldsymbol{e}}, {\boldsymbol{v}}) \right)+ O(\Delta x^2),\end{equation}

${\boldsymbol{e}}:=(1,0,0)$ is the unit vector in physical space, and similarly for the derivatives with respect to $y$ , $z$ , $v_x$ , $v_y$ and $v_z$ . The Vlasov equation discretized at spatial grid points using second-order accurate central differences can be expressed as

(2.5) \begin{align}\frac{\partial f({\boldsymbol{x}}, {\boldsymbol{v}}, t)}{\partial t} \approx& - \frac{{v}_x}{2\Delta x} \left( f(x+\Delta x, y, z, {\boldsymbol{v}}, t) - f(x-\Delta x, y, z, {\boldsymbol{v}}, t) \right) \nonumber \\[2pt]& - \frac{{v}_y}{2\Delta x} \left( f(x, y+\Delta x, z, {\boldsymbol{v}}, t) - f(x, y-\Delta x, z, {\boldsymbol{v}}, t) \right) \nonumber \\[2pt]& - \frac{{v}_z}{2\Delta x} \left( f(x, y, z+\Delta x, {\boldsymbol{v}}, t) - f(x, y, z-\Delta x, {\boldsymbol{v}}, t) \right) \nonumber \\[2pt]& - \frac{q}{m} \frac{\left( \textit{E}_x({\boldsymbol{x}}, t) + ({\boldsymbol{v}} \times {\boldsymbol{B}})_x({\boldsymbol{x}}, t) \right)}{2\Delta v}\left( \begin{array}{l} f({\boldsymbol{x}}, v_x+\Delta v, v_y, v_z, t) \\[2pt] - f({\boldsymbol{x}}, v_x-\Delta v, v_y, v_z, t) \end{array}\right) \nonumber \\[2pt]& - \frac{q}{m} \frac{\left( \textit{E}_y({\boldsymbol{x}}, t) + ({\boldsymbol{v}} \times {\boldsymbol{B}})_y({\boldsymbol{x}}, t) \right)}{2\Delta v}\left( \begin{array}{l} f({\boldsymbol{x}}, v_x, v_y+\Delta v, v_z, t) \\[2pt] - f({\boldsymbol{x}}, v_x, v_y-\Delta v, v_z, t) \end{array}\right) \nonumber \\[2pt]& - \frac{q}{m} \frac{\left( \textit{E}_z({\boldsymbol{x}}, t) + ({\boldsymbol{v}} \times {\boldsymbol{B}})_z({\boldsymbol{x}}, t) \right)}{2\Delta v}\left( \begin{array}{l} f({\boldsymbol{x}}, v_x, v_y, v_z+\Delta v, t) \\[2pt] - f({\boldsymbol{x}}, v_x, v_y, v_z-\Delta v, t) \end{array}\right).\end{align}

Introducing the vectorized value ${\boldsymbol{f}}(t)$ that represents the value of $f$ at each grid point, this equation can be transformed as follows:

(2.6) \begin{align} \frac {{\rm d} {\boldsymbol{f}}(t)}{{\rm d} t}= A(t) {\boldsymbol{f}}(t). \end{align}

Here ${\boldsymbol{f}}(t)$ is a $N^6$ -dimensional vector that is defined as

(2.7) \begin{align} {\boldsymbol{f}}_{j}(t) = f(j_x\Delta x,j_y\Delta x,j_z\Delta x,j_{v_x}\Delta v,j_{v_y}\Delta v,j_{v_z}\Delta v,t). \end{align}

We call a numerical scheme that discretizes only space using central differences, without discretizing time, a semidiscrete central difference scheme. When expressing it in terms of each spatial derivative term, $A(t)$ can be written as

(2.8) \begin{align} A(t) = A_x+A_y+A_z+A_{v_x}(t)+A_{v_y}(t)+A_{v_z}(t). \end{align}

Here $A_x$ is an $N^6 \times N^6$ square matrix, which can be expressed as the tensor product of each spatial matrix,

(2.9) \begin{align} A_x = -\frac {1}{2\Delta x}D_{\mathrm{per},x}\otimes I_y\otimes I_z\otimes \mathrm{diag}(v_x)\otimes I_{v_y}\otimes I_{v_z}, \end{align}

where $I_x$ is an $N \times N$ identity matrix and $\mathrm{diag}(v_x)$ is an $N \times N$ diagonal matrix with phase velocities as its diagonal elements. Here $D_{\mathrm{per},x}$ is an $N \times N$ central difference operator with the periodic boundary condition, which is defined as

(2.10) \begin{align} D_{\mathrm{per},x} = \begin{pmatrix} 0 & \quad 1 & \quad 0 & \quad \cdots & \quad 0 & \quad 0 & \quad -1 \\ -1 & \quad 0 & \quad 1 & \quad \cdots & \quad 0 & \quad 0 & \quad 0 \\ 0 & \quad -1 & \quad 0 & \quad \cdots & \quad 0 & \quad 0 & \quad 0 \\ \vdots & \quad \vdots & \quad \vdots & \quad \ddots & \quad \vdots & \quad \vdots & \quad \vdots \\ 0 & \quad 0 & \quad 0 & \quad \cdots & \quad 0 & \quad 1 & \quad 0 \\ 0 & \quad 0 & \quad 0 & \quad \cdots & \quad -1 & \quad 0 & \quad 1 \\ 1 & \quad 0 & \quad 0 & \quad \cdots & \quad 0 & \quad -1 & \quad 0 \end{pmatrix}\!.\nonumber\\[-12pt] \end{align}

From (2.10) and (2.9), it is clear that $A_x$ is an antisymmetric matrix. Here $A_{y}$ and $A_{z}$ are defined in the same manner and thus they are time-independent matrices.

On the other hand, $A_{v_x}$ is written as

(2.11) \begin{align} A_{v_x}(t) = &-\frac {1}{2\Delta v_x}\hat {E_x}(t)\otimes D_{\mathrm{per},v_x}\otimes I_{v_y}\otimes I_{v_z}\nonumber \\ &-\frac {1}{2\Delta v_x}\hat {B_z}(t)\otimes D_{\mathrm{per},v_x}\otimes \mathrm{diag}(v_y)\otimes I_{v_z}\nonumber \\ &+\frac {1}{2\Delta v_x}\hat {B_y}(t)\otimes D_{\mathrm{per},v_x} \otimes I_{v_y}\otimes \mathrm{diag}(v_z). \end{align}

Here $\hat {E_x}(t), \hat {B_y}(t)$ and $\hat {B_z}(t)$ are $N^3 \times N^3$ diagonal matrices with the electric and magnetic fields in the 3-D physical space as its diagonal elements. The differential operator in velocity space, $A_{v_x}, A_{v_y}$ , and $A_{v_z}$ are time-dependent. In this paper, however, we assume they are time-independent over a small time interval $[t, t + \Delta t]$ and define the effective Hamiltonian at each time step. That is, we use the time-independent Hamiltonian at the $n$ th time step as follows:

(2.12) \begin{align} \hat {H}_{t_n} := i A(t=t_n), \end{align}

where $t_n = n \Delta t$ and $n \in \{ 0, T/\Delta t - 1\}$ , with the simulation time $T$ . Because the matrix $A$ is antisymmetric, this effective Hamiltonian is a Hermitian.

Next, we define the effective quantum state as the sum of the normalized distribution function ${\boldsymbol{f}}_j$ encoded as the amplitude of the computational basis $|j\rangle$ corresponding to each grid point $j$ , normalized by ${\boldsymbol{f}}(t)\|$ ,

(2.13) \begin{align} |\psi (t)\rangle = \sum _{{\boldsymbol{f}}\in \mathcal{I}_6}\frac {{\boldsymbol{f}}_j(t)}{\|{\boldsymbol{f}}(t)\|}|j\rangle . \end{align}

In conclusion, we map the Vlasov equation to the Schrödinger equation with the time-dependent Hamiltonian,

(2.14) \begin{align} \frac {d }{{\rm d} t}|\psi (t)\rangle &= -i\hat {H}(t) |\psi (t)\rangle , \end{align}

and descretizing time, assuming that the Hamiltonian is constant over each time step (2.12). At each time step, the solution for the Schrödinger equation is given as

(2.15) \begin{align} |\psi (t_{n+1})\rangle = \mathrm{exp}\kern-0.5pt(-i\hat {H}_{t_n}\Delta t)|\psi (t_n)\rangle . \end{align}

Therefore, performing the time evolution iteratively, we finally obtain the general solution at any arbitrary time as

(2.16) \begin{align} |\psi (t)\rangle = \prod _{n=0}^{N_t-1}\mathrm{exp}\kern-0.5pt\big(-i\hat {H}_{t_n}\Delta t\big)|\psi (0)\rangle . \end{align}

Here $N_t$ is the number of time steps, and the step indices are defined as $n \in \{0, 1, 2, \ldots , N_t-1\}$ . This general solution is implemented using a quantum circuit.

2.2. Quantum singular value transformation

The QSVT is a quantum algorithm that applies polynomial transformations to the singular values of a given matrix (Gilyén et al. Reference Gilyén, Su, Low and Wiebe2019; Martyn et al. Reference Martyn, Rossi, Tan and Chuang2021). The QSVT is versatile and can be applied to various quantum algorithms, including Hamiltonian simulation. Specifically, it transforms the target function $\mathrm{exp}\kern-1pt(-i\hat {H}\Delta t)$ into an approximating polynomial $P_R$ of degree $2R+1$ , which is then implemented within a quantum circuit within the error tolerance $\epsilon$ .

In this study, we follow the discussions of Toyoizumi et al. (Reference Toyoizumi, Yamamoto and Hoshino2024) to formulate the necessary QSVT-based Hamiltonian simulation. First, block encoding is a method that encodes a matrix $\hat {H}$ as the upper-left block of a unitary matrix $U$ . For normalization factors $\eta , \alpha$ and the number of ancilla qubits $a$ , $U$ is called an $(\eta \alpha , a, 0)$ -block encoding of $\hat {H}$ if it satisfies the following condition:

(2.17) \begin{align} U = \begin{pmatrix} \frac {\hat {H}}{\eta \alpha } \quad \ \cdot \\ \cdot \qquad \cdot \end{pmatrix}\!, \end{align}

where the dot (.) denotes a matrix with arbitrary elements. Here, $\eta$ is defined as the normalization factor resulting from the Hadamard gates applied to the ancilla qubits, which depends on the dimensional degree of each problem set-up. For example, in the case of the 1D1V Vlasov equation, $\eta=\sqrt{2}^2$ and, the target matrix is embedded in the upper-left block of the unitary matrix, while the other blocks do not have physical significance.

When $\|\hat {H}\| \geqslant 1$ , it is necessary to use the normalization factor to normalize the matrix for block encoding to maintain unitarity. As discussed in Toyoizumi et al. (Reference Toyoizumi, Yamamoto and Hoshino2024), this implementation introduces a time parameter $\tau$ and adjusts the time step width to $\Delta t=\tau /\eta \alpha$ , thereby resolving any practical inconveniences in the implementation.

Next, using the implementation time parameter $\tau$ and the error $\epsilon$ , $U_{\mathrm{exp}}$ satisfies the following condition as a $(1, a+2, \epsilon )$ -block encoding of $\frac {1}{2}\mathrm{exp}\kern-2pt\left (-i ({(\hat {H})}/{(\eta \alpha )})\tau \right )$ :

(2.18) \begin{align} \left \Vert \frac {\mathrm{exp}\kern-2pt\left (-i({\hat {H}}/({\eta \alpha }))\tau \right )}{2}-\langle \text{ancilla}|_{0,1,2}\otimes U_{\mathrm{exp}} \otimes | \mathrm{ancilla} \rangle _{0,1,2}\right \Vert \leqslant \epsilon . \end{align}

Based on Toyoizumi et al. (Reference Toyoizumi, Yamamoto and Hoshino2024), figure 1 shows the $(1, a+2, \epsilon )$ -block encoding quantum circuit for $U_{\mathrm{exp}}$ . The phase factor $\phi \in \mathbb{R}^{d+1}$ in the figure corresponds to the degree $d$ of the approximation polynomial. The reason for the coefficient ${1}/{2}$ in $\mathrm{exp}\kern-1pt\big(-i ({(\hat {H})}/{(\eta \alpha )})\tau \big)$ is that, as shown in figure 1, we consider the attenuation of the amplitude of the solution state due to the superposition state of the ancilla qubit. According to Euler’s formula, we have

(2.19) \begin{align} \mathrm{exp}\kern-2pt\left (-ix\tau \right ) = \mathrm{cos}\kern-2pt\left (x\tau \right )-i\mathrm{sin}\kern-2pt \left (x\tau \right )\!, \end{align}

and we perform the Jacobi–Anger expansion as

(2.20) \begin{align} & \mathrm{cos}\kern-2pt \left (x \tau \right )=J_{0} (\tau )+2\sum _{k=1}^{\infty } (-1)^{k} J_{2k} (\tau )T_{2k} \left (x\right )\!,\nonumber \\ & \mathrm{sin}\kern-2pt \left (x \tau \right )=2\sum _{k=0}^{\infty } (-1)^{k} J_{2k+1} (\tau )T_{2k+1} \left (x \right )\!, \end{align}

where $J_{m} (\tau )$ is the $m$ th Bessel function of the first kind and $T_{k} (x)$ is the $k$ th Chebyshev polynomial of the first kind. We define the truncated series at an index $R$ by

(2.21) \begin{align} P_{R}^{\mathrm{cos}}(x) & :=J_{0} (\tau )+2\sum _{k=1}^{R} (-1)^{k} J_{2k} (\tau )T_{2k} (x),\nonumber \\ P_{R}^{\mathrm{sin} }(x) & :=2\sum _{k=0}^{R} (-1)^{k} J_{2k+1} (\tau )T_{2k+1} (x). \end{align}

Figure 1. The quantum circuit $U_{\mathrm{exp}}$ block encodes the time evolution operator $({1}/{2})\mathrm{exp}\kern-1pt(-i\hat {H}\Delta t)$ with $(1, a+2, \epsilon )$ -block encoding. Provided by Toyoizumi et al. Reference Toyoizumi, Yamamoto and Hoshino(2024).

According to Lemma 57 of Gilyén et al. (Reference Gilyén, Su, Low and Wiebe2019), the approximation polynomials $P_{R}^{\mathrm{cos}}(x)$ and $P_{R}^{\mathrm{sin} }(x)$ have an upper bound on the error under the following conditions:

(2.22) \begin{align} \left | \mathrm{cos}\kern-0.5pt (x\tau )-P_{R}^{\mathrm{cos}}(x)\right | & \leqslant \frac {5}{4}\left (\frac {e|\tau |}{4( R+1)}\right )^{( 2R+2)} ,\nonumber \\ \left | \mathrm{sin}\kern-0.5pt (x\tau )-P_{R}^{\mathrm{sin} }(x)\right | & \leqslant \frac {5}{4}\left (\frac {e|\tau |}{2( 2R+3)}\right )^{( 2R+3)} . \end{align}

Finally, substituting $\tau =\eta \alpha \Delta t$ and $x$ with ${(\hat {H})}/{(\eta \alpha )}$ , we get the $(1, a+2, \epsilon )$ -block encoding condition for $ ({1}/{2})\mathrm{exp}\kern-0.5pt\big(-i\hat {H}\Delta t\big)$ ,

(2.23) \begin{align} &\left \Vert \frac {\mathrm{exp}\kern-0.5pt(-i\hat {H}\Delta t)}{2}-\frac {\left ( P_{R}^{\mathrm{cos}}\left ( \hat {H}\Delta t\right ) -iP_{R}^{\mathrm{sin} }\left ( \hat {H}\Delta t\right )\right )}{2}\right \Vert \nonumber \\ &\quad{}\leqslant \frac {5}{8}\left (\frac {e|\eta \alpha \Delta t|}{4( R+1)}\right )^{( 2R+2)} +\frac {5}{8}\left (\frac {e|\eta \alpha \Delta t |}{2( 2R+3)}\right )^{( 2R+3)} ,\end{align}

(2.24) \begin{align} \left | \mathrm{cos}\kern-2pt \left ( \frac {\hat {H}}{\eta \alpha }\tau \right ) -P_{R}^{\mathrm{cos}}\left ( \frac {\hat {H}}{\eta \alpha }\tau \right )\right | \leqslant \varepsilon ,\left | \mathrm{sin}\kern-2pt \left ( \frac {\hat {H}}{\eta \alpha }\tau \right ) -P_{R}^{\mathrm{sin} }\left ( \frac {\hat {H}}{\eta \alpha }\tau \right )\right | \leqslant \varepsilon , \end{align}

where $0\lt \varepsilon \lt (1/e)$ and

(2.25) \begin{align} R( \tau ,\varepsilon ) & =\left \lfloor \frac {1}{2} r\left (\frac {e\tau }{2} ,\frac {5}{4} \varepsilon \right )\right \rfloor ,\nonumber \\ r( \tau ,\varepsilon ) & \leqslant \mathcal{O}( \tau +\log ( 1/\varepsilon )). \end{align}

Therefore, the number of queries for $U_{\mathrm{exp}}$ as $({1}/{2})\mathrm{exp}\kern-1pt\big(-i\hat {H}\Delta t\big)$ is given by

(2.26) \begin{align} R+R+1 & =2\left \lfloor \frac {1}{2} r\left (\frac {e\eta \alpha \Delta t}{2} ,\frac {5}{4} \varepsilon \right )\right \rfloor +1\\[-10pt]\nonumber \end{align}

(2.27) \begin{align} & \leqslant \mathcal{O}\left (\eta \alpha \Delta t +\log ( 1/\varepsilon )\right )\!.\\[12pt] \nonumber\end{align}

For more details see Gilyén et al. (Reference Gilyén, Su, Low and Wiebe2019). According to Toyoizumi et al. (Reference Toyoizumi, Yamamoto and Hoshino2024), amplitude amplification is necessary to cancel the factor $1/2$ in (2.23) and obtain $\mathrm{exp}\kern-1pt(-i\hat {H}\Delta t)$ . They present two types of QSVT-based algorithm to implement this: the oblivious amplitude amplification (Gilyén et al. Reference Gilyén, Su, Low and Wiebe2019) and the fixed-point amplitude amplification (Martyn et al. Reference Martyn, Rossi, Tan and Chuang2021). The total number of queries for $\mathrm{exp}\kern-1pt(-i\hat {H}\Delta t)$ , including amplitude amplification, is given by $3\left (2R+1\right )$ for oblivious amplitude amplification and $D\left (2R+1\right )$ for fixed-point amplitude amplification. The degree $D$ was given asymptotically in Gilyén et al. (Reference Gilyén, Su, Low and Wiebe2019) and Martyn et al. (Reference Martyn, Rossi, Tan and Chuang2021).

3.1. Quantum circuit for generating difference terms

The implementation of the quantum circuit is based on the coin operator and shift operator, which are the operational elements of the quantum walk. The coin operator provides transition probabilities, while the shift operator handles the vertex movement. See Douglas & Wang (Reference Douglas and Wang2009) for more details.

We define the rotation gate $R(\theta )$ along the $Y$ axis as

(3.1) \begin{align} R(\theta ) \equiv e^{-iY\arccos (\theta )} . \end{align}

In order to encode the coefficients of the difference term in (2.5) to the amplitudes of corresponding qubits, we apply the rotation gate to the $|0\rangle$ state as

(3.2) \begin{align} R(\theta )|0\rangle = \theta |0\rangle +\sqrt {1-|\theta |^2}|1\rangle , \end{align}

so that we obtain the desired state $|0\rangle$ with the coefficient $\theta$ encoded as its amplitude. In this encoding method, $\theta$ must be real and satisfy $|\theta | \leqslant 1$ . For general values including the cases with $|\theta | \gt 1$ , we need to normalize it using the factor $\alpha$ defined in (2.17). Specific choices of $\alpha$ will be represented in the following subsections.

Next, we perform the differentiation by applying $U_{\mathrm{Incr.}}$ and $U_{\mathrm{Decr.}}$ to the computational basis, which are defined as

(3.3) \begin{eqnarray} U_{\mathrm{Incr.}} | j \rangle = | j+1 \rangle ,\nonumber \\ U_{\mathrm{Decr.}} | j \rangle = | j-1 \rangle . \end{eqnarray}

These are unitary operators when we suppose the periodic boundary condition,

(3.4) \begin{eqnarray} U_{\mathrm{Incr.}} | N-1 \rangle = | 0 \rangle ,\nonumber \\ U_{\mathrm{Decr.}} | 0 \rangle = | N-1 \rangle , \end{eqnarray}

and can be implemented on a circuit as

(3.5)

Here, the solid and hollow dots indicate the 1 and 0 states of the control qubits.

As we correspond each computational basis with a grid point in our algorithm, incrementing or decrementing the computational basis corresponds to shifting the physical quantity at an arbitrary grid point by $\pm 1$ . In other words, with this operation, we obtain $f({\boldsymbol{x}}\pm {\boldsymbol{e}}\Delta x)$ .

3.2. Quantum circuit implementation for block encoding the effective Hamiltonian

In this subsection we enumerate specific effective Hamiltonians to use and circuits for the implementation.

Figure 2. Quantum circuit for the block encoding of the effective Hamiltonian of the 1-D Vlasov equation $\hat {H}_{1D}$ .

(i) Steady-state 1-D Vlasov equation (advection equation):

(3.6) \begin{equation} \frac {\partial f(x,t)}{\partial t} + v\frac {\partial f(x,t)}{\partial x} = 0. \end{equation}

When the phase velocity of the 1-D Vlasov equation is constant, it becomes equivalent to the steady-state advection equation. The advection equation describes the physics of the distribution function $f(x,t)$ propagating at the advection velocity $v$ while maintaining its waveform. By discretizing only the $x$ space using the central difference method to construct (2.6), we obtain $A = A_x$ . Note that normalization is performed using the normalization coefficient $\alpha$ to satisfy the condition $|\theta _v| \leqslant 1$ , where $\theta _v = v/2\Delta x \alpha = 1$ . Referring to Toyoizumi (Reference Toyoizumi2024), based on (2.9) and (2.12), the effective Hamiltonian can be expressed as

(3.7) \begin{align} \hat {H}_{1D} &= -i\theta _v D_{\mathrm{per},x}\nonumber \\ &= -i\theta _v \sum _{j_x=0}^{N-1}\left (|j_x\rangle \langle j_x+1|-|j_x\rangle \langle j_x-1|\right )\!. \end{align}

From figure 2, we can see that Hadamard gates are applied to the ancilla qubit, creating a superposition state. This means that the sum in (3.7) is constructed on the ancilla qubit’s $|0 \rangle _2$ state. Therefore, the amplitude is attenuated by a normalization factor $\eta =\sqrt {2}^2$ , corresponding to the number of Hadamard gates. Therefore, $(2\alpha ,2,0)$ th block encoding of the effective Hamiltonian is implemented in the quantum circuit as shown in figure 2.

(ii) The 1D1V Vlasov–Maxwell equations without a magnetic field:

(3.8) \begin{align} \frac {\partial f(x,v_x,t)}{\partial t} &+ v_x\frac {\partial f(x,v_x,t)}{\partial x} + \frac {q}{m}E(x,t)\frac {\partial f(x,v_x,t)}{\partial v_x} = 0,\nonumber\\ &\frac {\partial E(x,t)}{\partial t}=-\frac {q}{\epsilon _0}\int vf(x,v,t) \, {\rm d}v. \end{align}

The 1D1V Vlasov–Maxwell equations without a magnetic field describe the kinetic interaction between the plasma and the electric field. By discretizing only the phase space using the central difference method to construct (2.8), we obtain $A(t) = A_x + A_{v_x}(t)$ .

After discretizing the equation, we have the set of the coefficients $\{c_j\}$ with

(3.9) \begin{align} c_{j_{v_x}} &= \frac {v_{j_{v_x}}}{2\Delta x},\nonumber \\ c_{j_x} &= \frac {q}{m}\frac {E_x(x_{j_x},t)}{2\Delta v}. \end{align}

Therefore, we take the normalization factor $\alpha$ to satisfy the condition $|\theta | \leqslant 1$ as

(3.10) \begin{align} \alpha =\max \left (\left |\frac {v}{2\Delta x}\right |,\left |\frac {q}{m}\frac {E_x(x,t)}{2\Delta v}\right |\right )\!. \end{align}

Normalizing by this factor, we set $\theta$ as

(3.11) \begin{align} \theta _{j_{v_x}} &= \frac {v_{j_{v_x}}}{2\Delta x \alpha },\nonumber \\ \theta _{j_x} &= \frac {q}{m}\frac {E_x(x_{j_x},t)}{2\Delta v\alpha }. \end{align}

Based on (2.9), (2.11) and (2.12), the effective Hamiltonian can be described using $\theta$ as follows:

(3.12) \begin{align} \hat {H}_{1D1V} &= -\sum _{j_x=0}^{N-1}\sum _{j_{v_x}=0}^{N-1} \left ( \begin{aligned} &i\theta _{j_{v_x}}(|j_x\rangle \langle j_x+1|-|j_x\rangle \langle j_x-1|)\otimes |j_{v_x}\rangle \langle j_{v_x}|\\ &+i\theta _{j_x}|j_x\rangle \langle j_x|\otimes (|j_{v_x}\rangle \langle j_{v_x}+1|-|j_{v_x}\rangle \langle j_{v_x}-1|) \end{aligned} \right )\!. \end{align}

In a similar discussion as for the 1-D Vlasov case, the amplitude is attenuated by a normalization factor of $\eta =\sqrt {2}^4$ due to the four Hadamard gates. Therefore, the block encoding of the effective Hamiltonian for $(4\alpha ,3,0)$ is implemented as shown in figure 3. The electric field is updated at each time step by solving Ampere’s law from the Maxwell equations on a classical computer node. The straightforward implementations of $| \mathrm{phys} \rangle$ states and the controlled $R(\theta _{j_x})$ and $R(\theta _{j_{v_x}})$ gates require $O(N)$ gates, which does not provide a quantum advantage. According to Toyoizumi et al. (Reference Toyoizumi, Yamamoto and Hoshino2024), the number of gates can potentially be improved to $O(\text{poly}(\log (N)))$ assuming QRAM (Giovannetti, Lloyd & Maccone Reference Giovannetti, Lloyd and Maccone2008).

Figure 3. Quantum circuit for the block encoding of the effective Hamiltonian of the 1D1V Vlasov–Maxwell equations without magnetic field $\hat {H}_{1D1V}$ .

(iii) The 3D3V Vlasov–Maxwell equations with a magnetic field:

(3.13) \begin{align}\frac{\partial f({\boldsymbol{x}}, {\boldsymbol{v}}, t)}{\partial t}& + {\boldsymbol{v}} \boldsymbol{\cdot} \frac{\partial f({\boldsymbol{x}}, {\boldsymbol{v}}, t)}{\partial {\boldsymbol{x}}}+ \frac{q}{m} \left( {\boldsymbol{E}}({\boldsymbol{x}}, t) + {\boldsymbol{v}} \times {\boldsymbol{B}}({\boldsymbol{x}}, t) \right) \boldsymbol{\cdot} \frac{\partial f({\boldsymbol{x}}, {\boldsymbol{v}}, t)}{\partial {\boldsymbol{v}}}= 0,\end{align}

(3.14) \begin{align}\frac{\partial {\boldsymbol{E}}({\boldsymbol{x}}, t)}{\partial t}&= c^2 \nabla \times {\boldsymbol{B}}({\boldsymbol{x}}, t)- \frac{q}{\epsilon_0} \int {\boldsymbol{v}} f({\boldsymbol{x}}, {\boldsymbol{v}}, t)\, \mathrm{d}{\boldsymbol{v}}, \nonumber \\\frac{\partial {\boldsymbol{B}}({\boldsymbol{x}}, t)}{\partial t}&= - \nabla \times {\boldsymbol{E}}({\boldsymbol{x}}, t).\end{align}

The 3D3V Vlasov–Maxwell equations with a magnetic field fully reproduce the kinetic interaction between the plasma and the electromagnetic field. By discretizing only the 3D3V-phase space using the central difference method, we obtain (2.8). The normalization coefficient $\alpha$ is chosen to satisfy the condition $|\theta | \leqslant 1$ as

(3.15) \begin{align}\alpha = \max \left( \left| \frac{{\boldsymbol{v}}}{2\Delta x} \right|, \left| \frac{q}{m} \frac{{\boldsymbol{E}}({\boldsymbol{x}}, t) + {\boldsymbol{v}} \times {\boldsymbol{B}}({\boldsymbol{x}}, t)}{2\Delta v} \right|\right).\end{align}

This is defined as the sum of the maximum absolute values of the coefficients of each discretized term. Then, we set $\theta$ as

(3.16) \begin{align}\theta_{j_{v_x}} &= \frac{{v}_{j_x}}{2\Delta x\, \alpha}, \nonumber \\\theta_{j_{v_y}} &= \frac{{v}_{j_y}}{2\Delta x\, \alpha}, \nonumber \\\theta_{j_{v_z}} &= \frac{{v}_{j_z}}{2\Delta x\, \alpha}, \nonumber \\\theta_{{\boldsymbol{j}}_x, j_{v_y}, j_{v_z}} &= \frac{q}{m} \frac{\textit{E}_x({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)+ {v}_{y, j_{v_y}}\, \textit{B}_z({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)- {v}_{z, j_{v_z}}\, \textit{B}_y({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)}{2\Delta v\, \alpha}, \nonumber \\\theta_{{\boldsymbol{j}}_x, j_{v_x}, j_{v_z}} &= \frac{q}{m} \frac{\textit{E}_y({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)+ {v}_{z, j_{v_z}}\, \textit{B}_x({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)- {v}_{x, j_{v_x}}\, \textit{B}_z({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)}{2\Delta v\, \alpha}, \nonumber \\\theta_{{\boldsymbol{j}}_x, j_{v_x}, j_{v_y}} &= \frac{q}{m} \frac{\textit{E}_z({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)+ {v}_{x, j_{v_x}}\, \textit{B}_y({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)- {v}_{y, j_{v_y}}\, \textit{B}_x({\boldsymbol{x}}_{{\boldsymbol{j}}_x}, t)}{2\Delta v\, \alpha}. \nonumber\\[-18pt]\end{align}

Here, the index ${\boldsymbol{j}}_x:=(j_x,j_y,j_z)$ is defined as $\in \mathcal{I}_3:=[N]^{\otimes 3}$ . Therefore, based on (2.9), (2.11) and (2.12), the effective Hamiltonian can be expressed using $\theta$ as

(3.17) \begin{align}\hat{H}_{3D3V} = -\sum_{{\boldsymbol{j}} \in \mathcal{I}_6}^{N^6 - 1} \left(\begin{aligned} & i\theta_{j_{v_x}} \left( \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(x,+1)} - \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(x,-1)} \right) \\ + & i\theta_{j_{v_y}} \left( \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(y,+1)} - \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(y,-1)} \right) \\ + & i\theta_{j_{v_z}} \left( \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(z,+1)} - \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(z,-1)} \right) \\ + & i\theta_{{\boldsymbol{j}}_x,j_{v_y},j_{v_z}} \left( \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(v_x,+1)} - \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(v_x,-1)} \right) \\ + & i\theta_{{\boldsymbol{j}}_x,j_{v_x},j_{v_z}} \left( \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(v_y,+1)} - \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y-1,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(v_y,-1)} \right) \\ + & i\theta_{{\boldsymbol{j}}_x,j_{v_x},j_{v_y}} \left( \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(v_z,+1)} - \hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y-1,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(v_z,-1)} \right)\end{aligned}\right), \nonumber\\[-12pt]\end{align}

where we define

(3.18) \begin{align}\hat{{\boldsymbol{\textsf{T}}}}_{j_x,j_y,j_z,j_{v_x},j_{v_y},j_{v_z}}^{(x,\pm 1)}&:= |j_x\rangle \langle j_x \pm 1| \otimes |j_y\rangle \langle j_y| \otimes |j_z\rangle \langle j_z| \nonumber \\&\quad {} \otimes |j_{v_x}\rangle \langle j_{v_x}| \otimes |j_{v_y}\rangle \langle j_{v_y}| \otimes |j_{v_z}\rangle \langle j_{v_z}|,\end{align}

and similarly for the derivatives with respect to $y$ , $z$ , $v_x$ , $v_y$ and $v_z$ .

In a similar discussion as for the 1-D Vlasov case, the amplitude is attenuated by a normalization factor of $\eta =\sqrt {2}^8$ due to the eight Hadamard gates. Therefore, the $(16\alpha ,6,0)$ -block encoding of the effective Hamiltonian is implemented as shown in figure 4. The electric and magnetic fields are updated at each time step by solving Ampere’s law and Faraday’s law from the Maxwell equations on a classical computer. In the 3D3V case, we also assume that QRAM will be used to efficiently input the electromagnetic field information into the effective Hamiltonian at each time step. Finally, our proposed algorithm follows the flow shown in figure 5.

Figure 4. Quantum circuit for the block encoding of the effective Hamiltonian of the 3D3V Vlasov–Maxwell equations with magnetic field $\hat {H}_{3D3V}$ .

Figure 5. Our numerical simulation algorithm flow. The Vlasov equation is simulated on a quantum circuit using Hamiltonian simulation based on QSVT, while the velocity moment calculation and Maxwell equations are executed classically. Additionally, for comparison, we simulated by replacing the exponential matrix with Trotter decomposition of the quantum scheme and exact diagonalization on the classical node.

4.1. The 1-D advection test

To verify the properties of the QSVT-based Hamiltonian simulation as a numerical computation scheme, we conducted an advection test for the 1-D Vlasov equation with periodic boundary conditions. Additionally, for comparison, we implemented a first-order Trotter decomposition-based Hamiltonian simulation, referring to Sato et al. (Reference Sato, Kondo, Hamamura, Onodera and Yamamoto2024), and performed the exact diagonalization.

(i) Simulation conditions of 1-D advection test:

(4.1) \begin{align} f(x, t=0) = \begin{cases} 0 \quad \text{for } 0 \leqslant x \leqslant 2^{n_x-1}, \\ 1 \quad \text{otherwise}. \end{cases} \end{align}

From (2.13), we prepare the following state as an initial state:

(4.2) \begin{align} |\psi (t=0)\rangle = \frac {1}{\sqrt {2^{n_x-1}}}\sum ^{2^{n_x-1}-1}_{j=0}|j\rangle \otimes |1\rangle _{n_x-1}\otimes | \mathrm{ancilla} \rangle , \end{align}

here $n_x - 1$ qubits are used for preparing $|j\rangle$ state and four qubits are used for the ancilla state.

We set the number of $x$ qubits $n_x = 7$ (i.e. the number of grids $2^{n_x}=128$ ), the spatial interval $\Delta x = 1$ , the time range $0 \leqslant t \leqslant 18$ , the width of time step $\Delta t = 0.1$ , and the advection velocity $v = 1$ , resulting in the Courant number $\nu :=v\Delta t/\Delta x=0.1$ . The initial distribution function represents a rectangular wave. Other conditions include the degree of the QSVT approximation polynomial: $R = 14$ . The boundary condition is periodic.

Figure 6. The numerical results of each scheme for the 1-D Vlasov equation under the advection test conditions are shown. The solid line represents QSVT, the dashed line represents Trotter decomposition and the dotted line represents classical exact diagonalization. The red lines correspond to $t=0$ , the green lines to $t=9$ and the blue lines to $t=18$ , illustrating the time evolution.

Figure 7. The absolute errors between the numerical results of the quantum schemes and the classical exact diagonalization for the 1-D Vlasov equation under the advection test conditions are shown. The solid line represents the absolute error between QSVT and classical, and the dashed line represents the absolute error between Trotter decomposition and classical.

In this work, we execute quantum circuits on thestatevector_simulator, allowing us to obtain the values of $f_j(t)$ directly from the state vector. However, when real quantum devices are used, quantum circuits are executed multiple samples, and the values are obtained from the probability outcomes $p_j(t) = f_j^2(t) / \|f\|^2$ . Since $\|f\|$ is known initially and $f_j(t)$ is a positive real number, the desired value $f_j(t)$ can be straightforwardly extracted from $p_j(t)$ . The number of samples required to extract $f_j(t)$ within an error $\epsilon$ scales as $\mathcal{O}({(\|f\|^2)}/{\epsilon})$ .

Figure 6 shows the quantum computational results of the Hamiltonian simulations based on QSVT (solid line) and Trotter decomposition (dashed line), along with the numerical results of the classical exact diagonalization (dotted line). It shows the advection propagation with $v = 1$ at $t = 0$ (red), $t = 9$ (green) and $t = 18$ (blue). Numerical oscillations induced by the discretization using the central difference method are observed in the nonlinear gradient regions. These numerical oscillations can be mitigated by using the upwind difference method or by adding numerical viscosity. This can potentially be implemented by transforming the effective Schrödinger equation using the Schrödingerization method proposed by Hu et al. (Reference Hu, Jin, Liu and Zhang2024).

Next, figure 7 plots the absolute errors between QSVT and classical exact diagonalization (solid line), and between Trotter decomposition and the classical method (dashed line) (hereafter referred to as scheme error). The error of QSVT is seen to be proportional to $|f(x_j, t)|$ and time $t$ at a given time $t$ at grid points $x_j$ . On the other hand, the error in Trotter decomposition is influenced by the high-frequency components of the numerical oscillations. This indicates that significant differences in numerical results can arise depending on which quantum scheme is used. We want to emphasize that, as we move towards practical applications, it is essential to select the appropriate quantum scheme based on the physical problem settings. Moreover, since the results of Trotter decomposition and classical exact diagonalization are consistent with those presented in Sato et al. (Reference Sato, Kondo, Hamamura, Onodera and Yamamoto2024), the validity of our code is assured.

4.2. The 1D1V two-stream instability test

Next, we focused on the two-stream instability test. Historically, in the development of classical Vlasov solvers, the two-stream instability is one of the most fundamental nonlinear problems for the 1D1V Vlasov–Maxwell equations of plasma. Whether this phenomenon can be accurately reproduced under appropriate conditions serves as a criterion for evaluating the numerical computation scheme (Filbet & Sonnendrücker Reference Filbet and Sonnendrücker2003; Crouseilles, Mehrenberger & Sonnendrücker Reference Crouseilles, Mehrenberger and Sonnendrücker2010; Qiu & Shu Reference Qiu and Shu2011).

We utilized a single node with an A100 GPU 40 GB from the Computing Infrastructure Center at the University of Tsukuba, executing the quantum circuit simulator with Qiskit-Aer-GPU. The solutions to the Maxwell equations, computed on a classical node, were sequentially encoded into the effective Hamiltonian of the QSVT-based Hamiltonian simulation to solve the Vlasov–Maxwell system. Thus, it should be noted that this algorithm is a quantum–classical hybrid algorithm.

(i) Simulation conditions of 1D1V two-stream instability based on Qiu & Shu (Reference Qiu and Shu2011):

(4.3) \begin{align} f(x,v,t=0) &= \frac {2}{7\sqrt {2\pi }} \left (1 + \beta \frac {\mathrm{cos}\kern-0.5pt (2kx) + \mathrm{cos}\kern-0.5pt (3kx)}{1.2} + \mathrm{cos}\kern-0.5pt (kx)\right ) \mathrm{exp}\kern-1.5pt\left (-\frac {v^{2}}{2}\right )\!, \nonumber \\ E(x,t=0) &= 0. \end{align}

From (2.13), we prepare the following state as an initial state:

(4.4) \begin{align} |\psi (t=0)\rangle = \sum ^{2^{n_x-1}-1}_{j_x=0}\sum ^{2^{n_v-1}-1}_{j_v=0}\frac {f(j_x\Delta x,j_v\Delta v,t=0)}{\|f(t=0)\|}|j_v\rangle _v \otimes |j_x\rangle _x \otimes | \mathrm{ancilla} \rangle _{0,1,2,3,4,5}. \end{align}

Here six qubits are used for the ancilla state. We set the number of qubits for the $x$ and $v$ axis as $n_x = n_v = 5$ , (i.e. simulation domain is defined the 1D1V phase space $32\times 32$ ), the time range $0 \leqslant t \leqslant 2500$ , the time step $\Delta t = 1.89$ , the $x$ space resolution $\Delta x = 0.39$ , the $V$ space resolution $\Delta v = 0.31$ , the initial perturbation amplitude parameter $\beta = 0.01$ and the initial wavenumber parameter $k = 0.5$ . The grid size was chosen to achieve approximately $ \nu \approx 1$ . The initial distribution function represents relative velocity in phase space. Other parameters were set as follows: the degree of the QSVT approximation polynomial $R = 9$ , charge $q = 1$ , mass $m = 1$ and permittivity $\epsilon _0 = 1$ . The boundary condition for the $x$ physical space is periodic, while the boundary condition for the $v$ velocity space is fixed. The implementation method for fixed boundary conditions in QSVT is described in Appendix A. For simplicity, the current density is assumed to be the first moment of the distribution function for a single species of ions equation (3.8).

Figure 8. The quantum computational numerical results of the QSVT-based Hamiltonian simulation scheme for the 1D1V Vlasov–Maxwell equations under the two-stream instability conditions. The distribution function in the $x$ – $v$ phase space is shown at real time $T = 53$ , with the vertical axis representing the velocity space and the horizontal axis representing the physical space.

Figure 8 shows the results of the two-stream instability at $T = 53$ for the 1D1V Vlasov–Maxwell equations using QSVT. The occurrence of phase mixing is evident. Figure 9 shows the results obtained using classical exact diagonalization, and the vortex formation in the phase space is almost identical. However, focusing on the colour bar, the values of the distribution function in QSVT are generally smaller than those in the classical method. This discrepancy arises from the quantum scheme error (approximation error) in QSVT’s $ \mathrm{exp}(-i\hat {H}\Delta t)$ under the condition of (2.23), leading to a difference from the true values. Nevertheless, as shown in figure 7, the quantum scheme error in QSVT is proportional to the magnitude of the distribution function, thus not significantly affecting the vortex formation. Therefore, the results between QSVT’s (figure 8) and the classical method’s (figure 9) appear almost identical.

These results are also reasonably consistent with the vortex size and phase mixing depiction in the numerical results of the two-stream instability using the classical high-precision Vlasov scheme MPP SL DG method by Qiu & Shu (Reference Qiu and Shu2011). However, our results show the onset of numerical oscillations around the $x$ intervals $[-4,-2]$ and $[2,4]$ . Over long-term evolution, these numerical oscillations are likely to induce different physical phenomena. As discussed in the 1-D advection test, to ensure monotonicity, it is necessary to either introduce upwind differencing or incorporate numerical viscosity to suppress numerical oscillations and improve accuracy for nonlinear developments.

Figure 9. The numerical results of the classical exact diagonalization scheme for the 1D1V Vlasov–Maxwell equations under the two-stream instability conditions. The distribution function in the $x$ – $v$ phase space is shown at real time $T = 53$ , with the vertical axis representing the velocity space and the horizontal axis representing the physical space.

Figure 10. The norms of distribution function and electric field of the QSVT Hamiltonian simulation scheme for the 1D1V Vlasov–Maxwell equations under the two-stream instability conditions. Panels (a) and (b) show the L1 norm and the L2 norm, respectively.

Figure 11. The norms of distribution function and electric field of the classical exact diagonalization scheme for the 1D1V Vlasov–Maxwell equations under the two-stream instability conditions. Panels (a) and (b) show the L1 norm and the L2 norm, respectively.

Figure 10 shows the time evolution of the L1 and L2 norms of the distribution function and the electric field obtained from the QSVT Hamiltonian simulation. Due to the operator norm error in each QSVT step, the L1 and L2 norms of the distribution function decrease exponentially (as does the probability). However, this decrease can be mitigated by increasing the degree of the QSVT approximation polynomial. In other words, the error can be made arbitrarily small by increasing computational cost. For more details, see § 5. Figure 11 shows the time evolution of the L1 and L2 norms of the distribution function and the electric field obtained from classical exact diagonalization. Since the Vlasov equation preserves the L1 and L2 norms of the distribution function in the case of the periodic boundary condition in this work, unlike the QSVT results, the classical exact diagonalization without operator norm error preserves these norms.