Phase-space grid for a boundary-value problem

What it fills

• Discretize a PDE into one linear system

  Replace every continuous variable at once — space together with time, or space together with velocity — so that the whole problem becomes a single matrix equation. Nothing is left to march: the grid values at every recorded point are unknowns of one system, solved in one go.

When it applies

Outgoing, non-reflecting boundary conditions are imposed at both spatial edges to avoid artifacts from reflected waves, and an artificial diffusivity $\eta$ is imposed in velocity space so that the grid resolution can be kept low. Novikau, Dodin and Startsev are explicit that the diffusivity is not free: it complicates the block-encoding and increases the condition number of the resulting matrix. Their scheme is validated against an analytical solution for homogeneous plasma only, and they describe the system as a minimal problem chosen as a testbed rather than a practical application.

Requires

Every step this method names moves its route along, so there is nothing it needs alongside them.

Example

given  linearized Vlasov-Ampere perturbation g(x,v) and electric field E(x)
       for a 1D, collisionless, neutral Maxwellian electron plasma,
       driven by a source current j^(S)(x)                             (Eq. 5)
       spatial grid  x_j = j*h,        h = xmax/qx,   j = 0..Nx-1      (Eq. 15)
       velocity grid v_k = -vmax+k*dv, dv = 2*vmax/qv, k = 0..Nv-1     (Eq. 16)
       Nx = 2^nx, Nv = 2^nv,   qx = Nx-1, qv = Nv-1                    (Eq. 17, Sec. III)

requires  outgoing (non-reflecting) boundary conditions at both spatial
    edges, so that waves reflected from the box do not pollute the
    solution -- this is what the factor zeta_bc below encodes  (Sec. II B, Eq. 14)

optional  an artificial diffusivity eta >= 0 in velocity space. It is a free
    parameter, NOT a precondition: the paper imposes it "to avoid numerical
    artifacts and keep the grid resolution reasonably low", and its own
    comparison against analytical theory is run at eta = 0 (nx=9, nv=8).
    Turning it on is a trade -- the paper says it complicates the later
    block-encoding and somewhat raises the condition number of A
                                        (Sec. II B before Eq. 12; Sec. V, Fig. 2)

# --- reduce to a boundary-value problem: no time axis remains --------------
g, E  ~  exp(-i*omega0*t)                                               (Eq. 8)
    # assumption: a single real drive frequency, not derived from the
    #   dynamics -- this is what removes d/dt and turns the initial-value
    #   Vlasov-Ampere system into a boundary-value one, so nothing here
    #   is marched in time                                              (Sec. II B)
i*omega0*g - zeta_bc*v*(dg/dx) + eta*(d2g/dv2) - v*H*E = 0               (Eq. 13a)
i*omega0*E + integral( v*g, dv ) = j^(S)                                 (Eq. 13b)
    # H = F/T: the background Maxwellian F over its temperature T, both
    #   frozen coefficients of the linearization, not unknowns   (Sec. II B, after Eq. 9b)
    # zeta_bc = 0 for an incoming wave at that edge, 1 otherwise -- this
    #   is what makes the boundaries outgoing / non-reflecting          (Eq. 14)
    # with eta = 0 this is Eq. (9a) unmodified; Eq. (12) is (9a) plus the
    #   eta*d2g/dv2 term                                                (Eq. 9a, 12)

# --- central differences in the bulk ----------------------------------------
dg/dx   |_{j,k} = sigma*(g_{j+1,k} - g_{j-1,k}),        sigma = 1/(2h)   (Eq. 19a)  j = 1..qx-1
d2g/dv2 |_{j,k} = beta*(g_{j,k+1} -2*g_{j,k}+ g_{j,k-1}), beta = 1/dv^2  (Eq. 20a)  k = 1..qv-1

# --- one-sided edges from a Lagrange interpolating polynomial --------------
dg/dx   |_{0,k}  = sigma*(-3*g_{0,k}   +4*g_{1,k}    -g_{2,k})           (Eq. 19b)
    # 2nd-order Lagrange polynomial, spatial edges                      (Sec. III)
dg/dx   |_{qx,k} = sigma*(3*g_{qx,k}   -4*g_{qx-1,k} +g_{qx-2,k})        (Eq. 19c)
d2g/dv2 |_{j,0}  = beta*(2*g_{j,0}     -5*g_{j,1}+4*g_{j,2}-g_{j,3})     (Eq. 20b)
    # 3rd-order Lagrange polynomial, velocity edges                     (Sec. III)
d2g/dv2 |_{j,qv} = beta*(2*g_{j,qv}    -5*g_{j,qv-1}+4*g_{j,qv-2}-g_{j,qv-3})  (Eq. 20c)

# --- fold into one coefficient per grid point, then the two field eqns -----
P_{j,k} = i*omega0 + zeta_bc_{j,k}*(delta_{j,0}-delta_{j,qx})*3*v_k*sigma
              - p_sign_k * 2*eta*beta                                   (Eq. 22)
    # (delta_{j,0}-delta_{j,qx}) and p_sign_k = 1-2*(delta_{k,0}+delta_{k,qv})
    #   are bookkeeping for the sign flips the one-sided stencils above
    #   carry at the spatial and velocity edges -- not new physics      (Eq. 22, 26)
    # at eta = 0 the last term vanishes and the construction is unchanged
    #   otherwise; nothing downstream depends on eta being nonzero      (Eq. 22)

P_{j,k}*g_{j,k} + sum_i P_x_{j,k,i}*g_{i,k} + sum_i g_{j,i}*P_v_{i,k}
    - v_k*H_{j,k}*E_j = 0                                               (Eq. 21)
    # P_x, P_v package the stencils above as coefficients of neighboring
    #   x- and v-grid points respectively                               (Eq. 25, 26)
i*omega0*E_j + sum_k  v_k*g_{j,k}  = j^(S)_j                             (Eq. 27)
    # the dv factor from the continuous integral (13b) is already gone
    #   here -- it was absorbed once, up front, by renormalizing g and
    #   F themselves (dv*g -> g, dv*F -> F)                             (Eq. 18)

# --- stack g_{j,k} and E_j into one vector, read off A and b ---------------
psi = stack( g_{j,k} for all j,k ;  E_j at k=0, zero elsewhere )        (Eq. 28, 29)
A = [[ F,   C^E ],
     [ C^f, S   ]]                                                      (Eq. 32)
    # F holds the P_{j,k} coefficients and the edge/bulk stencils above,
    #   as three block types F^L, F^B, F^R for left edge / bulk / right
    #   edge spatial points                                             (Eq. 33-36)
    # C^E, C^f are block-diagonal: the v*H*E and Ampere sum-over-v
    #   couplings between g and E respectively                          (Eq. 37, 38)
    # S = i*omega0 * I, diagonal                                        (Eq. 39)
b = stack( 0 for every g_{j,k} slot ; j^(S)_j at k=0, zero elsewhere )   (Eq. 30, 31)

return  the linear system  A*psi = b                                    (Eq. 1)
# no time-stepping loop exists in this method -- the fixed-frequency
#   assumption above already converted "march forward in time" into
#   "solve one linear system", which is the point of casting it as
#   Eq. (1) in the first place                                          (Sec. I, II B)
# discretization error itself is not bounded anywhere in the paper;
#   only the condition number's sensitivity to eta is quantified, and
#   only at one resolution (nx=7, nv=5, omega0=1.2, eta=0.002 against
#   eta=0) -- not a general cost bound                                  (Sec. V)
# validated only against an analytical solution for HOMOGENEOUS plasma;
#   the paper calls the whole setup a minimal testbed, not yet a
#   practical target                                                    (Sec. V, Sec. I)

Cost, as the source states it

The paper states the price of the diffusivity as measured condition numbers at one resolution rather than as a general bound: at $n_x = 7$, $n_v = 5$ and $\omega = 1.2$, the matrix has $\kappa_A = 8.844 \times 10^4$ with $\eta = 0.002$ against $\kappa_A = 3.489 \times 10^4$ without it. No bound on the discretization error itself is stated, and none is recorded here.

Implementations

• KIN1D1D, the authors' own CUDA assembly and solve of $\mathbf{A}\boldsymbol{\psi} = \boldsymbol{b}$

  • Encoding of linear kinetic plasma problems in quantum circuits via data compression

  The paper cites its own solver. Reference 23 of arXiv:2403.11989 reads "Classical modeling of the electrostatic kinetic plasma problem., https://github.com/QuCF/KIN1D1D (2024)", and it is cited exactly once in the body, in Section V, on the sentence "The numerical calculations were performed by inverting the matrix (32) using the sparse-QR-factorization-based method provided in CUDA toolkit cuSOLVER". The problem it exists to answer is the one that section opens with — "To test our discretization scheme, we have performed classical simulations for homogeneous plasma ($n = T = 1$), which facilitates comparison with the analytic theory described in Appendix A" — so this is not a demonstration written up after the paper but the instrument the discretization scheme was checked with, and it produces the matrix $\mathbf{A}$ of $\mathbf{A}\boldsymbol{\psi} = \boldsymbol{b}$ rather than a propagator, because there is no time axis left to march. Every commit but the repository-creation stub is Ivan Novikau's, the paper's first author — the oldest, `e9ecbca9` ("Initial commit", 2023-04-28), is authored by the `QuCF` account itself — and the newest is "SVD" on 2024-04-27, six weeks after the preprint. Two things a reader should carry in before opening it. The repository is shared with a neighbouring problem and its three self-descriptions do not agree: the README says "Modeling of 1D-1D stationary electrostatic kinetic plasma system", GitHub's repository description says "Anomalous skin effect modeling using 1D real and 1D velocity spaces", and the class header comment in `include/kin1d.cuh` says "Linear kinetic waves in a 1-D electron plasma (skin effect)"; what is recorded here is the path `src/main.cu` actually takes, which is the electrostatic one. And the committed configuration is one of this paper's runs rather than all of them — `diff = -0.002` is the $\eta = 0.002$ setting of Figs. 4 and 5, not Fig. 2's $\eta = 0$ — so it is the shipped outputs, not the source alone, that tie the artefact to each figure.

  `kin_outgoing_form1` in `src/main.cu` constructs a `KWout_diff` — declared in `include/kin_outgoing_diff.cuh`, deriving from `KW` in `include/kin1d.cuh` — and then calls `init_device`, `set_background_profiles`, `save_xv_background_profs_to_hdf5`, `form_rhs`, `save_rhs_to_hdf5`, `form_sparse_matrix`, `save_matrices`, `solve_system`, `recheck_result` and `save_result`, in that order and once each, with `dd.svd()` sitting commented out between `form_sparse_matrix` and `save_matrices`; there is no loop over time anywhere in the driver, which is the shape of the method rather than an omission. The `KW` constructor lays the joint grid down exactly as Eqs. (15) and (16): `dd_.Nx = 1 << nx` and `dd_.Nv = 1 << nv`, then `dd_.h = dd_.xmax / (dd_.Nx - 1)` and `dd_.dv = 2.*dd_.vmax / (dd_.Nv - 1)`, with `dd_.A.N = dd_.Nx * dd_.Nv * dd_.Nvars` and `dd_.Nvars = 2` for the $g$ and $E$ halves. The stencils are two CUDA kernels in `device/dev_outgoing_diff.cuh`, each launched one thread per velocity index and one block per spatial index — `Diff_set_sparse_matrix_rows` over `dd_.Nx+1` blocks of `dd_.Nv` threads writes the CSR row pointers, `Diff_set_matrix_values` over `dd_.Nx` blocks of `dd_.Nv` threads writes the values. Three device helpers carry the velocity-space second derivative $\eta\partial_v^2$ — not the sign of $v$, which is decided elsewhere — and all three scale by `idv = dev_dd_.diff/(dev_dd_.dv*dev_dd_.dv)`: `bulk_vel` writes the triple $(-1, +2, -1)$ at columns `ir - 1`, `ir`, `ir + 1`, which is the central second difference of Eq. (20a); `left_vel_boundary` writes $(-2, +5, -4, +1)$ at `ir` through `ir + 3` and `right_vel_boundary` its mirror at `ir - 3` through `ir`, which are Eqs. (20b) and (20c), the third-order Lagrange one-sided pair, under the convention that the paper's $\eta$ is the negative of the code's `diff`. The spatial derivative rides on `ih = v1/(2.*dev_dd_.h)`, the paper's $v_k\sigma$: bulk spatial points take $+v_k\sigma$ at `ir - Nv` and $-v_k\sigma$ at `ir + Nv`, and at the two edges the three-point one-sided coefficients of Eqs. (19b) and (19c) appear as `ih3` on the diagonal with `-ih4` and `ih` beyond it on the left, and `-ih`, `ih4`, `-ih3` on the right. $\zeta^{\rm bc}$ is not stored as a coefficient at all: the `if(iv < Nvh)` branches simply omit the advective entries for the incoming half of velocity space at each edge, so Eq. (14) is implemented as absent matrix entries rather than as a multiplying zero — which is also why the incoming rows are shorter, and why the row-offset kernel `Diff_set_sparse_matrix_rows` carries the separate per-row strides `n_out = Nd + 1 + Ndv` and `n_in = 2 + Ndv`. `KWout_diff::compute_Nnz`, one file over, never sees those two and gets the total from a single closed form, `3 * Nv * (Nx - 2) + 4 * Nv + (Ndv*(Nv-2) + 2*Ndv_bo) * Nx` plus `3 * Nv * Nx`, with `Ndv_bo = 3` and `Ndv = 2` set in the `KWout_diff` constructor. The three remaining blocks are written inline at the end of `Diff_set_matrix_values`: `-v1*FB[ir]` for $\mathbf{C}^{E}$, a bare `v1` for $\mathbf{C}^{f}$, and an imaginary `w` for the diagonal $\mathbf{S}$. `FB` is filled in `device/dev_profiles_rhs.cuh` as `coef * exp(-v2/(2*T1))` after `coef = coef * dev_dd_.dv / T1`, its own comment reading "dv * Maxwellian / T" — the paper's $H = F/T$ with the $\Delta v\,F \to F$ half of Eq. (18) already folded in, which is why $\mathbf{C}^{E}$'s entry needs no $\Delta v$ beside it; the other half of Eq. (18), $\Delta v\,g \to g$, is separately why $\mathbf{C}^{f}$'s Ampère entry is the bare `v1`, and `FB` appears in no other matrix block. `KW::solve_system` at `include/kin1d.cuh` line 381 calls the member `KW::solve_sparse_system` at line 776, which at line 782 calls `LA::solve_sparse_system` — `cusolverSpZcsrlsvqr`, double-precision complex sparse QR, reached through the `ycuInverse` macro at `include/mix.h` line 55 — and wraps it in a `YTimer` whose result is printed and never written to a file. `KWout_diff::svd`, the branch `main.cu` comments out, instead hands a densified copy to `LA::cond_number_cusolver`, built on `cusolverDnZgesvd`.

  No dataset is read; both the background and the right-hand side are generated. `set_background_profiles` with `id_profile = "flat"` sets `T_[ii] = 1.0` and `den_[ii] = 1.0` at every spatial point, which is the homogeneous $n = T = 1$ case Section V validates against, and `build_background_distribution` fills the Maxwellian from them; `tanh2` and `exp` profiles exist in the same function and are not the path taken. The source is `init_rhs_gauss` in `device/dev_profiles_rhs.cuh`, writing `b[ir].y = coef_ch * coef_norm * exp(-dr2/ds2)` with `coef_ch = dev_dd_.w` and `coef_norm = 1.0` into the $E$ slots only and into the imaginary part only, i.e. $j^{(S)}_j = \mathrm{i}\omega_0\exp[-(x_j-x_0)^2/2\Delta_S^2]$. Only $n_x$, $n_v$ and $\omega_0$ are read from `argv`; the rest of the run is hard-coded in `src/main.cu` — `Lx = 100`, `Lv = 4`, `source_x0 = 50.`, `source_ds = 1.` at lines 41 to 44, `id_prof = "flat"` at line 55, and above them `Tref = 1.0e4 * cc.ev_` and `den_ref = 1e12`. Those match Section V's stated $x_{\max} = 100$, $n_x = 9$, $v_{\max} = 4$, $n_v = 8$, source at $x_0 = 50$ with $\Delta_S = 1.0$ and homogeneous plasma; `diff = -0.002` on line 45 does not, being the $\eta = 0.002$ of Figs. 4 and 5, and the zero-diffusivity runs were made by editing this file, as the live `folder_to_save` ending `diff0002` and its commented-out `diff0` alternative on the next line show. The outputs are what ships: `simulations/EVM/classical-sims/` in the companion repository github.com/QuCF/QuCF holds 31 HDF5 files running from `out_3_3_w1.2_Lx100_Lv4_flat.hdf5` to `out_9_8_w1.2_Lx100_Lv4_flat.hdf5` — one at $n_x = n_v = 3$ and then $n_x = 4$ through $9$ against $n_v = 4$ through $8$ — beside `w08/` with three files at $\omega_0 = 0.8$ and `w12-zero-diff/` with 25 at $\eta = 0$. Read on 2026-08-27, `out_7_5_w1.2_Lx100_Lv4_flat.hdf5` records `basic/launch-path` as `/media/work/docs/codes/QuCF/KIN1D1D/build` and `basic/date-of-simulation` as `08-10-2023 13:00:19`, `parameters` as `w` $= 1.2$, `diff` $= -0.002$, `source_x0` $= 50.0$, `source_dx` $= 1.0$, and grids of $128$ and $32$ points spanning $[0, 100]$ and $[-4, 4]$; and it carries the whole assembled system rather than only the answer — `matrices/A-values`, `A-columns`, `A-rows`, `A-N`, `A-Nnz`, `profiles/F`, `profiles/b` and `result/psi`. The two files matching Fig. 2's stated parameters are here too, under different names and off different builds: `w12-zero-diff/out_9_8_w1.2_Lx100_Lv4_flat.hdf5` has `parameters/diff` of $-0.0$ at `matrices/A-N` $= 262144$, while `w08/out_9_8_w0.8_Lx100_Lv4_flat.hdf5` carries no `parameters/diff` and no `matrices` group at all, so it came off the plain `KWout` path rather than the `KWout_diff` one recorded here.

  github.com/QuCF/KIN1D1D, branch `main`, CUDA C++ with GitHub reporting the language as Cuda, no tags and no releases, and no forks or stars. There is no licence, and that is a checked absence rather than an assumption: the repository root holds only `.gitignore`, `README.md`, `build/`, `device/`, `include/` and `src/`, GitHub's licence endpoint returns 404 for the repository, and no file under `src/`, `include/` or `device/` carries a copyright or licence header — a case-insensitive search for "copyright", "licence", "license" and "GPL" across the nine files read here returned nothing — so there is no header text to quote and a reader has no grant to reuse it. Line counts read 2026-08-27 at HEAD, commit `b3c157c5` ("SVD", Ivan Novikau, 2024-04-27, the newest commit in the repository): `include/kin1d.cuh` 919 lines, `include/mix.h` 885, `device/dev_outgoing_diff.cuh` 382, `include/LA.cuh` 260, `include/kin_outgoing.cuh` 139, `include/kin_outgoing_diff.cuh` 100, `src/main.cu` 87, `device/dev_profiles_rhs.cuh` 56, `include/kinetic_data.h` 51. `device/dev_outgoing_diff.cuh` has two commits, `b3c157c5` and `f0997cc2` ("new", 2023-10-05); `include/kin1d.cuh` has four, those two plus `504ceb15` ("EVM", 2024-01-28) and `5c2b08a0` ("initial configuration", 2023-05-05). The four files to open are `device/dev_outgoing_diff.cuh` — `Diff_set_sparse_matrix_rows` at line 14, `left_vel_boundary` at 84, `bulk_vel` at 118, `right_vel_boundary` at 145, `Diff_set_matrix_values` at 184 — then `include/kin1d.cuh` for `class KW` at line 10, the grid at lines 111 and 112 and `solve_sparse_system` at 776, `include/kin_outgoing_diff.cuh` for `class KWout_diff` at line 6, `form_sparse_matrix` at 29 and `compute_Nnz` at 72, and `src/main.cu` for `kin_outgoing_form1` at line 24. The shipped outputs live in a second, likewise unlicensed repository, github.com/QuCF/QuCF, branch `main`, HEAD `6c6f957b` of 2025-01-11, under `simulations/EVM/`; the post-processing lives in a third, github.com/QuCF/scripts-py, under `jupyter-notebooks/EVM/`.

  The paper reports its validation as a figure and one sentence rather than as an error norm, and the sentence is not about the runs the rest of this entry follows. Fig. 2 is the high-resolution, zero-diffusivity comparison: it plots $\mathrm{Re}\,E$ and $\mathrm{Im}\,E$ "computed numerically (blue) and analytically (red) using Eq. (A17)" for $\omega_0 = 1.2$ and $\omega_0 = 0.8$, and its caption ends "In both cases, $n_x = 9$, $n_v = 8$, and $\eta = 0$" — with no error figure attached to it. The paper's one worded agreement claim belongs to a lower-resolution, diffusive pair instead: "Artifacts become noticeable at lower resolution (Figs. 4 and 5) but can be suppressed by introducing artificial diffusivity $\eta$ in the velocity space [Eq. (12)]. Such simulations are demonstrated in Figs. 4 and 5 for $n_x = 7$, $n_v = 5$, and $\eta = 0.002$. As seen in Fig. 4, their results are in good agreement with the analytical solution." Those are the parameters the committed `src/main.cu` carries and the parameters of the file read below, so the conditioning quoted here and the agreement claimed there are about the same run. The only scalars drawn from any of these simulations are conditioning — "if one takes $n_x = 7$, $n_v = 5$, $\eta = 0.002$, $\omega = 1.2$, the condition number of the resulting matrix is $\kappa_A = 8.844 \times 10^4$ (i.e. $\log_2 \kappa_A = 16.4$). Without the diffusivity and with the same resolution, the condition number is $\kappa_A = 3.489 \times 10^4$ (i.e. $\log_2 \kappa_A = 15.1$)" — and no discretization-error bound and no run time is stated anywhere. Both numbers survive being recomputed from the shipped artefacts, which is worth saying because the matrix itself is in the files: reading `matrices/A-*` out of `out_7_5_w1.2_Lx100_Lv4_flat.hdf5` and its `w12-zero-diff/` twin on 2026-08-27 gives an $8192 \times 8192$ matrix with $32960$ nonzeros in both cases, and `scipy.sparse.linalg.svds` at $k = 1$ returns $\sigma_{\max} = 16.366$ while a shift-invert `eigsh` on $\mathbf{A}^{H}\mathbf{A}$ returns $\sigma_{\min} = 1.8505 \times 10^{-4}$ at `parameters/diff` of $-0.002$, hence $\kappa_A = 8.844 \times 10^4$, against $\sigma_{\min} = 4.6906 \times 10^{-4}$ and $\kappa_A = 3.489 \times 10^4$ at `parameters/diff` of $-0.0$; the shipped `result/psi` satisfies its own system to $\max\lvert \mathbf{A}\boldsymbol{\psi} - \boldsymbol{b}\rvert = 5.0 \times 10^{-16}$ and $7.2 \times 10^{-16}$ respectively. The same figures are hard-coded as a scan table in `jupyter-notebooks/EVM/condition_number.ipynb` of github.com/QuCF/scripts-py, whose $n_v = 5$ rows read `sp_norm` `1.637e+01`, `min_sing` `1.851e-04` and `cond_nu` `8.844e+04` at $n_x = 7$ and `eta = -0.002`, against `min_sing` `4.691e-04` and `cond_nu` `3.489e+04` at `eta = 0.0`. Nothing quantum runs in this artefact: the solve is `cusolverSpZcsrlsvqr` on whatever CUDA device the binary is launched on, the wall time the solver measures is printed to stdout and saved in no committed file, and neither the repository nor the paper records a GPU model.

• Classiq's `vlasov_ampere` notebook, the $(x,v)$ stencils rebuilt twice

  • Implementation of a quantum linear solver for the Vlasov-Ampere equation
  • Encoding of linear kinetic plasma problems in quantum circuits via data compression

  `applications/plasma/vlasov_ampere` in Classiq's open library is the same discretization derived again by people who did not write it. The `description` field of the directory's `vlasov_ampere.metadata.json` opens "This demonstration was created in collaboration with the authors of the paper on encoding linear kinetic plasma problems in quantum circuits." and goes on for three more sentences about QSVT and Qmod; the notebook's own first cell says it "is based on the paper" of arXiv:2403.11989 "and was created in collaboration with its authors"; and the paper documenting it, arXiv:2507.22257, closes its acknowledgments with "We thank Ilya Dodin and Ivan Novikau for helpful discussions" and states its scope as "We implement a quantum linear solver for the one-dimensional Vlasov-Ampere equation, following the model presented in Novikau et. al." What earns it a place under this method rather than under the block-encoding one beside it is that the notebook rebuilds the grid and the stencils from scratch, twice: once as a plain NumPy matrix builder used as the classical reference, and once as Qmod functions whose structure deliberately mirrors it — the paper's own contrast being that "the original work by Novikau et al. also presented the block-encoding construction, albeit using low-level details and techniques. In contrast, the Qmod approach provides a more transparent formulation, highlighting the connection between the classical matrix and its quantum data representation." It is a reduced version of the system, and the reduction is silent: the artificial velocity diffusivity $\eta$ and with it the entire $\partial_v^2$ operator of Eqs. (20) are absent from the notebook — neither of the two governing equations it states carries the term, and `get_block_encoding` returns a three-term sum with no fourth — so what is implemented throughout is the $\eta = 0$ case and the third-order Lagrange velocity-edge stencils never appear.

  The classical builder is the notebook's second-to-last code cell, cell index 70. `get_advective_mat(nx, nv, dx, v_max, cyclic=False)` forms the $x$ derivative as `np.diag(np.ones(2**nx - 1), k=1) - np.diag(np.ones(2**nx - 1), k=-1)` and then overwrites the two edge rows with `boundary = np.pad([-3, 4, -1], (0, dx_mat.shape[1] - 3))` and `dx_mat[-1, :] = -np.flip(boundary)`, which is Eqs. (19b) and (19c) written out; it takes `np.kron(dx_mat / (2 * dx), np.diag(v_amplitudes))` for $v_k\sigma$ and then applies $\zeta^{\rm bc}$ as a row mask rather than a factor, `xi = 1 - np.kron(x_values == 0, v_amplitudes > 0) - np.kron(x_values == max_x, v_amplitudes <= 0)` followed by `advective[xi == 0] = 0`. `get_off_diag_mat` fills the couplings, and the two are not the same shape: `A[i*dim_v : (i+1)*dim_v, i*dim_v] = vec_col` writes `v_H_amplitudes` down the first column of each diagonal block, `B[i*dim_v, i*dim_v : (i+1)*dim_v] = vec_row` writes `v_amplitudes * dv` across the first row of each, and they are assembled as `top = np.hstack((np.zeros_like(A), -A))` and `bottom = np.hstack((B, np.zeros_like(B)))` — the paper's $\mathbf{C}^{E}$ of Eq. (37) with its minus sign in place, a column, and $\mathbf{C}^{f}$ of Eq. (38), a row. `get_block_encoding` returns `1j * w0 * np.eye(mat_size) - np.kron(np.diag([1, 0]), advective) + off_diag`, and `solve_problem` closes with `np.linalg.solve(be, source)`. The quantum side is the same matrix as Qmod `qfunc`s over a `BEData` struct of `v: QNum[N_V, SIGNED, 0]`, `x: QNum[N_X]` and a single `E: QBit` that switches between the $g$ and $E$ halves. `derivative_dirichlet_be` is the two-diagonal bulk part as an `lcu([1, -1], ...)` over `inplace_add(-1, extended_qnum)` and `inplace_add(1, extended_qnum)`, the $x$ register bound together with a `flag` qubit so overflow is caught rather than wrapped; `derivative_boundary_min_be` prepares the single row `BC_VALUES = 0.5 * np.array([-3, 3, -1, 0])`, normalised in place by its own `BC_NORM`, with `inplace_prepare_amplitudes`, and `derivative_boundaries_be` reverses the row order and adds a $\pi$ phase under control of the $x$ register's top qubit, so one prepared row serves both edges; `derivative_be` combines the two with `lcu([1, BC_NORM], ...)`, and their sum, at the scale `DERIVATIVE_TERM_FACTOR = 2 * (1 + BC_NORM)` restores, is the $(-3, 4, -1)$ and $(1, -4, 3)$ edge rows the markdown above it draws out as an explicit matrix identity. `v_be` is the diagonal $v$ factor through `assign_amplitude_table(lookup_table(lambda n: n / (2 ** (v.size - 1)), v), v, ind)`; `zeta_be` is two comparisons, `flag ^= (x == 0) & (v > 0)` and `flag ^= (x == (2**x.size - 1)) & (v <= 0)`, with a comment noting the two conditions are mutually exclusive so no disjunction is needed; `advective_be` multiplies the three and kills the lower-right block with `be.block.block_flag ^= be.data.E`. The two off-diagonal terms differ in the quantum build the way they differ in the classical one, and the order of the two statements is what carries it: `force_term_be` fires `flag ^= v != 0` first and then loads `-v_H_amplitudes`, emptying columns, while `current_term_be` applies `invert(lambda: load_v_vector(v))` — the adjoint of loading `v_amplitudes` — and fires `flag ^= v != 0` after it, emptying rows, which is the "first row in each block" shape $\mathbf{C}^{f}$ has in Eq. (38). `off_diag_be` selects between them on the `E` qubit and rebalances the two with `equalize_amplitude`, and `full_be` is `lcu(LCU_COEFFS, [advective_be, off_diag_be, IDENTITY], ...)` with `LCU_COEFFS = [-ADVECTIVE_TERM_FACTOR, OFF_DIAG_FACTOR, 1j * W_0, 0]` — the minus on the advective term and the $\mathrm{i}\omega_0$ on the identity being Eq. (13a) read straight off, and the trailing zero a pad to the two-qubit `lcu` register.

  No dataset; every grid and every coefficient is computed in the notebook. The parameters are one cell — `Temperature = 1`, `N = 1`, `DS = 3`, `X_0 = 50`, `W_0 = 0.8`, `X_MIN = 0`, `X_MAX = 100`, `V_MIN = -4`, `V_MAX = 4`, `N_X = 3`, `N_V = 3` — an $8 \times 8$ phase-space grid against the paper's $2^9 \times 2^8$, so $\Delta x = 100/7 \approx 14.3$ Debye lengths, which is why the notebook says of its own answer that "the problem solution is not physical, though we show it for educational purposes". Two departures from the paper's grid are worth naming because nothing in the notebook flags either. The velocity grid is `np.roll(np.linspace(-1, 1 - 2 ** (-N_V + 1), 2**N_V) * V_MAX, 2**N_V // 2)`, which at `N_V = 3` evaluates to $[0, 1, 2, 3, -4, -3, -2, -1]$: a two's-complement ordering for the signed `QNum`, spanning $[-v_{\max},\ v_{\max} - 2v_{\max}/2^{n_v}]$ with step $2v_{\max}/2^{n_v}$, not the paper's symmetric $v_k = -v_{\max} + k\Delta v$ with $\Delta v = 2v_{\max}/(N_v-1)$ of Eq. (16). And the `dv` that `get_block_encoding` computes and passes into `get_off_diag_mat`, where it multiplies the Ampère row, is the paper's `(2 * v_max) / (2**nv - 1)`, equal to $8/7$ at `N_V = 3` — the same value the parameter cell's `DV` carries into `CURRENT_TERM_FACTOR` on the quantum side — while the spacing of the `v_amplitudes` grid it weights is $1$. The source is built twice and the two builds are not equal. The classical reference's `get_source_term` writes `j_amplitudes = 1j * w0 * np.exp(-((x_coordinates - x0) ** 2) / (2 * ds**2))` into the $E$ slots only, the same shape KIN1D1D generates and at the same $x_0 = 50$, but with $\Delta_S = 3$ against the paper's $1.0$. The quantum half's `prepare_source_term` instead loads `j_amplitudes = W_0 * np.exp(-((x_coordinates - X_0) ** 2) / (2 * DS**2))`, with the trailing comment `# ignore the i phase`, divided by its own `SOURCE_NORM_FACTOR` — so a notebook whose point is to build one system twice differs between its two halves, deliberately, by the factor $\mathrm{i}$ and by a normalisation at exactly this point, while the markdown cell above the quantum one still writes the source with its $\mathrm{i}\omega_0$.

  github.com/Classiq/classiq-library, branch `main`, MIT licensed — `LICENSE.txt` opens "MIT License" and "Copyright (c) 2024 Classiq Technologies" — 2031 stars, latest tag `1.26.0`, all read 2026-08-27. The artefact is the directory `applications/plasma/vlasov_ampere/`, holding `vlasov_ampere.ipynb`, `vlasov_ampere.qmod`, `vlasov_ampere_qiskit.ipynb`, `vlasov_ampere.metadata.json`, `vlasov_ampere.synthesis_options.json` and a PDF of the accompanying paper. Same date: `vlasov_ampere.ipynb` is 402058 bytes, nbformat 4.9, 73 cells of which 31 are code, saved from a Python 3.11.4 kernel and carrying its execution outputs; `vlasov_ampere.qmod` is 9645 bytes of synthesised Qmod holding the `BELcu`, `BE`, `BEData` and `BEBlock` structs and the expanded `qfunc`s below them; `vlasov_ampere_qiskit.ipynb` is 54800 bytes and 37 cells, pinned to `qiskit==2.1.1`, and says of itself that it "implements the same model that appears under" `vlasov_ampere.ipynb`, "using qiskit" and is "the code used to collect the benchmarking data in the paper". The notebook has 23 commits; it was added by `8f703e97` ("with plasma notebook, without execution", Or Samimi Golan, 2025-07-31), its last substantive change is `ef893cd3` of 2026-07-07, which migrated the notebook to Classiq's new execution API, and HEAD for the file is `28c4e7a8` ("Title-case all markdown headings across 136 notebooks", 2026-07-30). The functions a reader should open are `get_advective_mat`, `get_off_diag_mat`, `get_block_encoding`, `get_source_term` and `solve_problem` in cell 70 for the classical matrix, and `derivative_dirichlet_be`, `derivative_boundary_min_be`, `derivative_boundaries_be`, `derivative_be`, `v_be`, `zeta_be`, `advective_be`, `force_term_be`, `current_term_be`, `off_diag_be` and `full_be` for the Qmod one.

  The discretization reproduces exactly, and that is the one result in this directory that belongs to this method rather than to the block encoding above it. Running the notebook's own `get_advective_mat` and `get_block_encoding` here on 2026-08-27 at `N_X = N_V = 3` gives Eq. (19b) to the digit: with $\Delta x = 100/7$, the $x = 0$, $v = -4$ row carries $0.42$, $-0.56$, $0.14$ at $x$-columns $0$, $1$, $2$, which is $(-3, 4, -1)\,v_k/(2\Delta x)$, and its $x = 7$, $v = 3$ mirror carries $0.105$, $-0.42$, $0.315$, which is Eq. (19c); the $\zeta^{\rm bc}$ mask then empties the $x = 0$, $v > 0$ rows and the $x = 2^{n_x}-1$, $v \le 0$ rows entirely, leaving zero nonzeros in each. The $\le$ is the code's, not the paper's: `get_advective_mat` masks on `v_amplitudes <= 0` at the right edge and `zeta_be` on `v <= 0`, where Eq. (14) and the notebook's own markdown both write $v < 0$. The extra row that empties is $v_k = 0$, whose advective entries carry a factor $v_k$ and vanish either way, so at this grid the two conditions give the same matrix. Everything else the notebook stores is resource counting on a simulator. Synthesising the block-encoding model with its reference register prints "Num qubits: 25" and the full QSVT model prints "Num qubits: 19"; the QSVT phase fit at `DEGREE = 180` prints `min singular value: 0.01856162144286599, max_singular value: 0.8282275541075902` and `Max relative error value: 0.02046509592876211`, these being singular values of the extracted, normalised block-encoded matrix rather than of $\mathbf{A}$ itself; and the execution is `calculate_state_vector(qprog_qsvt, filters={"block": 0})`, a statevector simulation, with the recovered $\lvert E(x)\rvert$ plotted in cell 68 and the `np.linalg.solve` reference plotted in cell 71 with the builder cell between them. No hardware run appears anywhere in the directory, and the qiskit sibling names no backend either. The accompanying paper's numbers are resource counts too: "The construction of the block encoding using the high-level Qmod language involves $n_x + n_v + 1$ logical qubits for the data variable and additional 8 logical qubits for the block variable", the benchmarked single QSVT step "operates on $n_x + n_v + 10$ logical qubits", and "Classiq approach can reduce the CX-counts by two orders of magnitudes, compared to the baseline Qiskit result" across grid sizes from $(2^3, 2^3)$ to $(2^6, 2^6)$ in its Fig. 1. Its discussion concedes the limit this method's own record already carries: "The specific problem considered here is not expected to exhibit quantum advantage, primarily due to the scaling of the matrix condition number with problem size."

What it needs

Nothing below this — it bottoms out here.

Other ways to fill the same slot

Different approaches

• Forward-time, centre-space (FTCS)

  Take the forward difference in time and the three-point central difference in space, then stack the resulting one-step relations for every timestep into a single block lower-bidiagonal system whose unknowns are the grid values at all recorded times together.

In the Atlas

• Encoding of linear kinetic plasma problems in quantum circuits via data compression

  Encode a linear kinetic plasma problem — modeling electrostatic linear waves, driven by a spatially localized external current, in a one-dimensional Maxwellian electron plasma — into a quantum circuit that solves the resulting linear system.

Sources

• Encoding of linear kinetic plasma problems in quantum circuits via data compression↗Ivan Novikau, Ilya Y. Dodin, Edward A. Startsev · 2024theory yes · simulation ? · hardware ?