# Quantized tensor networks for solving the Vlasov-Maxwell equations

Erika Ye<EMAIL_ADDRESS>Applied Mathematics and Computational Research
Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Rd., Berkeley, CA
94720  Nuno F. G. Loureiro<EMAIL_ADDRESS>Plasma Science and Fusion
Center, Massachusetts Institute of Technology, 77 Massachusetts Av.,
Cambridge, MA 02139

###### Abstract

While the Vlasov-Maxwell equations provide an ab-initio description of
collisionless plasmas, solving them is often impractical due to high
computational costs. In this work, we implement a semi-implicit Vlasov-Maxwell
solver utilizing the quantized tensor network (QTN) framework. This framework
allows one to efficiently represent and manipulate low-rank approximations of
high-dimensional data sets. As a result, the cost of the solver scales
polynomially with parameter $D$ (the so-called bond dimension), which is
directly related to the error associated with the low-rank approximation. By
increasing $D$, convergence to the dynamics that the solver would obtain
without any low-rank approximation is guaranteed. We find that for the 2D3V
test problems considered here, a modest $D=64$ appears to be sufficient for
capturing the expected physics, despite the simulations using a total of
$2^{36}$ grid points and thus requiring $D=2^{18}$ for exact calculations.
Additionally, we utilize a QTN time evolution scheme based on the Dirac-
Frenkel variational principle, which allows us to use larger time steps than
that prescribed by the Courant-Friedrichs-Lewy (CFL) constraint. As such, the
QTN format appears to be a promising means of approximately solving the
Vlasov-Maxwell equations with significantly reduced cost.

##  I Introduction

“Quantum-inspired” algorithms loosely describe a set of classical algorithms
that utilize some concepts from their quantum counterparts. They often obtain
comparable speed-ups so long as some error with respect to the true solution
can be tolerated and the original problem of interest has some (potentially
unknown) underlying structure that can be capitalized upon. In the context of
solving partial differential equations (PDEs), quantum-inspired algorithms
utilize a quantum-like representation of the classical dataset. One can then
use classical methods for quantum simulation to approximately solve the PDE.
One natural choice is the tensor network framework, a popular tool in the
quantum chemistry and condensed matter communities Schollwöck (2011); Vidal
(2003); White (1993); Chan and Head-Gordon (2002); White (2005) since it may
allow for efficient representation and manipulation of the quantum-like
dataset, so long as the dataset is low rank.

This concept was first introduced several years ago, and is referred to as
quantics or quantized tensor trains (QTTs) because the data is artificially
folded into a high-dimensional tensor representing many small “quants” of data
Khoromskij (2011); Dolgov  _et al._ (2012); Oseledets (2009, 2010). (Tensor
trains are tensor networks with a 1-D geometry.) It was also proven that
several one-dimensional functions (exponential, trigonometric, polynomial) and
common operators (e.g., differential operators) on a uniform grid can be
efficiently represented in the QTT format Khoromskij (2011); Kazeev and
Khoromskij (2012); Kazeev  _et al._ (2013a). It follows that one can often
obtain a loose bound on the QTT rank for other smooth functions by expanding
them in terms of exponential or polynomial basis sets Dolgov  _et al._ (2012).
However, non-smooth functions can have efficient representations too; the
Heaviside step function can be represented with just rank two, triangle
functions with rank three Gourianov  _et al._ (2022). The underlying intuition
is that functions with some degree of scale locality, meaning that features on
long length scales are relatively independent of behavior at small length
scales, can be represented as a tensor network with low rank. Since it is
believed that turbulent systems also exhibit some degree of scale locality,
tensor networks may even be able to efficiently capture turbulent dynamics.

There have been several works that solve the Navier-Stokes equations using the
quantized tensor train (QTT) format Gourianov  _et al._ (2022); Kiffner and
Jaksch (2023); Kornev  _et al._ (2023); Peddinti  _et al._ (2023) with
reasonable success. In this work, we are interested in solving the Vlasov
equation, a (6+1)-dimensional PDE that provides an ab initio description of
plasma dynamics in the presence of self-consistent electromagnetic fields
computed via coupling to Maxwell’s equations. Being able to perform these
calculations would greatly advance our understanding of many astrophysical
phenomena and fusion energy systems. While the coupled Vlasov-Maxwell system
is the gold standard for plasma simulation, it is seldom used due to
prohibitively high computational costs. A commonly used alternative are
Lagrangian or Monte-Carlo solvers (particle-in-cell), but those face
challenges in managing statistical noise and are limited in the range of
physical scales they can cover. There are a few works that investigate the use
of low-rank solvers (e.g., Refs. Kormann, 2015; Einkemmer and Lubich, 2018;
Allmann-Rahn  _et al._ , 2022; Einkemmer, 2023); they show that, often, only
modest rank is required to perform common test problems with sufficient
accuracy. However, while algorithmically similar, those investigations are
fundamentally different from the QTN methods investigated here in that they do
not make low-rank approximations within each dimension.

Prior work on solving the Vlasov equation with QTTs assumed electrostatic
fields, and only considered 1D1V (1 spatial axis, 1 velocity axis)
calculations Ye and Loureiro (2022). Here, we consider higher dimensional
plasmas in 2D3V (though the algorithm can be used for 3D3V calculations) in
the presence of electromagnetic fields. We first start with a brief overview
of the Vlasov equation and QTTs, and then introduce a quantized tensor network
(QTN) geometry that we expect to be a more efficient representation for high-
dimensional data than the standard tensor train geometry. In the following
section, we detail our semi-implicit finite-difference Vlasov solver in the
QTN format. Most notably, we discuss how one can modify the problem such that
positivity of the distribution function is preserved despite the low-rank
approximation, and introduce a QTN time-evolution scheme that allows one to
take larger time steps than allowed by traditional grid-based methods.
Finally, we close with numerical simulations of the Orszag-Tang vortex and the
GEM magnetic reconnection problem in 2D3V and a discussion on our
observations.

##  II Problem Overview

The collisionless Vlasov equation is given by

$\frac{\partial f_{s}}{\partial
t}+\textbf{v}_{s}\cdot\nabla_{\textbf{x}}f_{s}+\textbf{F}_{s}\cdot\nabla_{\textbf{v}}f_{s}=0,$
(1)

where $f$ is the distribution function with coordinates in real space (x) and
velocity space (v). The subscript $s$ denotes the particle species (e.g., ions
or electrons) that is being described. F denotes the force that the particles
are subject to. In the case of plasma dynamics, the force of interest is the
Lorentz force,

$\textbf{F}_{s}=q_{s}\left(\frac{\textbf{v}_{s}\times\textbf{B}}{m_{s}}+\textbf{E}\right),$
(2)

where $q_{s}$ is the charge and $m_{s}$ is the mass of particle species $s$, B
is the magnetic field, and E is the electric field.

For an electromagnetic system, one propagates the electric and magnetic fields
forward in time via Maxwell’s equations

$\displaystyle\frac{\partial\textbf{B}}{\partial t}$
$\displaystyle=-\nabla\times\textbf{E},$ (3)
$\displaystyle\frac{\partial\textbf{E}}{\partial t}$
$\displaystyle=c^{2}\nabla\times\textbf{B}-\frac{1}{\varepsilon_{0}}\textbf{J},$
(4)

where $c$ is the speed of light, $\varepsilon_{0}$ is the vacuum permittivity,
and
$\textbf{J}=\sum_{s}q_{s}\int\,d\textbf{v}\,\textbf{v}f_{s}(\textbf{x},\textbf{v})$
is the current density. The current density is the source of the nonlinearity,
since the state of the plasma affects the electromagnetic force that is
exerted on itself.

Eqs. (1-4) are the equations we wish to solve in our Vlasov-Maxwell solver.

##  III Quantized Tensor Networks

Suppose a 1-dimensional function $f(x)$ is discretized into $N=d^{L}$ grid
points. In the quantized tensor train format Oseledets (2009); Khoromskij
(2011); Dolgov (2014), the data is reshaped into an $L$-dimensional tensor
with size $d$ along each dimension, and encoded into a 1-D chain of $L$
tensors:

$\displaystyle f(x_{i})$ $\displaystyle\cong f(i_{1},i_{2},...,i_{L})=$
$\displaystyle\sum_{\alpha_{1}=1}^{r_{1}}...\sum_{\alpha_{L-1}=1}^{r_{L-1}}M^{(1)}_{\alpha_{1}}(i_{1})M^{(2)}_{\alpha_{1},\alpha_{2}}(i_{2})\ldots
M^{(L)}_{\alpha_{L-1}}(i_{L}),$ (5)

where $M^{(p)}$ is the $p^{\text{th}}$ tensor in the chain with one “physical”
bond ($i_{p}$) of size $d$ that is used to access the actual data, and two
“virtual” bonds ($\alpha_{p-1},\alpha_{p}$) of size $r_{p-1}$ and $r_{p}$ that
connect the tensor to its two neighbors. (The exceptions are $M^{(1)}$, which
only has one virtual bond ($\alpha_{1}$) of size $r_{1}$, and $M^{(L)}$, which
also only has one virtual bond ($\alpha_{L-1}$) of size $r_{L-1}$.) Low-rank
approximations in the QTT format limit the size of the virtual bonds to a
maximum value of $D$, which is often referred to as the bond dimension of the
TT. (The terms bond dimension and rank will be used interchangeably.) If the
original dataset can be represented as a QTT of small bond dimension (ideally
$D\sim\log(N)$) with only some small tolerable amount of error, the dataset is
described as compressible and can be efficiently represented as a QTT.

The indices $i_{1},i_{2},...i_{L}$ encode the grid points in $x$. For example,
if $d=2$, one possible encoding is a binary mapping, such that the
$i^{\text{th}}$ grid point is indexed by its binary string
$(i_{1},i_{2},...,i_{L})_{2}$. With this mapping, each tensor sequentially
correspond to finer grid resolutions. (For concrete examples of common
functions in the QTT format, we refer the reader to page 56 of Ref. Dolgov,
2014.) Another possible encoding is a signed register, in which one takes the
two’s complement of the physical indices conditioned on the first index, such
that $i\cong(i_{1},i_{1}+i_{2},...,i_{1}+i_{L})_{2}$. This mirror-like mapping
may be more efficient for symmetric functions Ripoll (2021).

The QTT representation can be considered “quantum-inspired” since mapping
classical data onto a quantum computer could be done in a similar fashion
Ripoll (2021); Lubasch  _et al._ (2020); Jaksch  _et al._ (2023). In this
case, each tensor would correspond to the state of one qubit.

Operators acting on a function $f(x)$ can also be expressed in the QTT format,
resulting in the QTT-Operator

$\displaystyle O(x^{\prime},x)$ $\displaystyle\cong
O(o_{1},o_{2},...,o_{L};,i_{1},i_{2},...,i_{L})=$
$\displaystyle\sum_{\alpha_{1}=1}^{r_{1}}...\sum_{\alpha_{L-1}=1}^{r_{L-1}}M^{(1)}_{\alpha_{1}}(o_{1},i_{1})M^{(2)}_{\alpha_{1},\alpha_{2}}(o_{2},i_{2})\ldots$
$\displaystyle\hskip 85.35826pt\ldots M^{(L)}_{\alpha_{L-1}}(o_{L},i_{L})$ (6)

where each tensor $M^{(p)}$ now has two physical bonds ($o_{p},i_{p}$) that
are used to access the output $x^{\prime}$ and input $x$ of the discretized
operator $O$.

In the case of $K$-dimensional data sets with $K>1$, one can consider several
different ansatzë. Recent works utilizing the QTN format thus far have only
used the 1-D tensor network (tensor train) geometry. One can opt to append
tensors corresponding to different dimensions sequentially (Fig. 1(a)) Ye and
Loureiro (2022); Kiffner and Jaksch (2023); Kornev  _et al._ (2023); Peddinti
_et al._ (2023) or in an interleaved fashion (Fig. 1(b)) Gourianov  _et al._
(2022). The sequential ordering is perhaps more intuitive, but interactions
between different dimensions can be long-range. For example, any correlation
between $x$ and $z$ are long-range and must be transmitted through all tensors
corresponding to the $y$ grid points. With the interleaved geometry, in which
tensors are ordered by grid scale, interactions between different dimensions
are typically not as long-range. However, this ordering is less efficient for
separable functions and operators.

(a) sequential, serial$x$$y$$z$

(b) sequential, interleaved$\hat{x}$$\hat{y}$$\hat{z}$

(c) comb TTN$\hat{x}$$\hat{y}$$\hat{z}$

Figure 1: Tensor network diagrams for three-dimensional data. (a) Tensor train
with tensors for different dimensions appended sequentially. (b) Tensor train
with tensors for different dimensions combined in an interleaved fashion. (c)
Comb-like tree tensor network with tensors for different dimensions on
separate branches.

This work uses a comb-like tree tensor network ansatz (Fig. 1(c)) to represent
a $K$-dimensional dataset,

$\displaystyle f(x_{1,i},x_{2,i},\dots,x_{K,i})\cong$ $\displaystyle
f(i^{(1)}_{1},\dots,i^{(1)}_{L},i^{(2)}_{1},\dots,i^{(2)}_{L},\dots,i^{(K)}_{1},\dots,i^{(K)}_{L})=$
$\displaystyle\hskip
14.22636pt\sum_{\gamma_{1}=1}^{r_{1}}...\sum_{\gamma_{K-1}=1}^{r_{K-1}}B^{(1)}_{\gamma_{1}}(i_{1},...i_{L})B^{(2)}_{\gamma_{1},\gamma_{2}}(j_{1},...,j_{L})\ldots$
$\displaystyle\hskip 85.35826pt\ldots
B^{(K)}_{\gamma_{L-1}}(i^{(K)}_{1},...,i^{(K)})\,,$ (7)

where $B^{(k)}$ is a tensor train for the $k^{\text{th}}$ dimension

$\displaystyle B^{(k)}_{\gamma_{k},\gamma_{k+1}}(i_{1},...,i_{L})=$
$\displaystyle\hskip
14.22636pt\sum_{\beta_{k}}\tilde{M}^{(k)}_{\gamma_{k},\gamma_{k+1},\beta_{k}}\biggl{(}\sum_{\alpha_{k,1}}...\sum_{\alpha_{k,L-1}}M^{(k,1)}_{\alpha_{k},\alpha_{(k,1)}}(i^{(k)}_{1})$
$\displaystyle\hskip
56.9055ptM^{(k,2)}_{\alpha_{(k,1)},\alpha_{(k,2)}}(i^{(k)}_{2})\ldots
M^{(k,L)}_{\alpha_{(k,L-1)}}(i^{(k)}_{L})\biggr{)}\,.$ (8)

In this geometry, the grid points along each dimension are represented by 1-D
tensor trains (the branches) and each TT is connected to each other via
another tensor train (the spine; darker green). In the above equations, tensor
$\tilde{M}$ is a spine tensor and does not include a physical index. To
simplify notation, we assumed that all dimensions have the same number of grid
points. However, this is not a necessary requirement. This comb-like tensor
network (QTC) has also been referred to as the QTT-Tucker format Dolgov and
Khoromskij (2013). For simplicity, this work uses the same bond dimension for
all virtual (internal) bonds, though it may be advantageous to use different
bond dimension for the bonds in the branches ($\alpha_{(k,i)}$), along the
spine ($\gamma_{k}$), and between the spine to the branches ($\beta_{k}$).

The motivation for this ansatz is to reduce the distance along the tensor
network between different dimensions, while maintaining an efficient
representation for product states. Additionally, all algorithms for tensor
trains can be extended to the comb geometry in a straightforward fashion,
since the comb geometry contains no loops and one can still define a center of
orthogonality. The computational costs of relevant algorithms are summarized
in Table 1. The scaling of most algorithms are more expensive for the comb
geometry because of the spine tensors, though in practice, this may be offset
by increased representability, allowing one to use smaller bond dimension.

Procedure | QTT | QTC-Spine
---|---|---
Advection in x |  |
| QTT-Operator and QTT-State multiplication
---
| $\mathcal{O}(D^{3}D_{W}^{2}d+D^{2}D_{W}^{3}d^{2})$
---
\+ ($Dd\times Dd$) eigenvalue problem
|  $\mathcal{O}(S^{4}S_{W})$ to $\mathcal{O}(S^{4}S_{W}^{4})$
---
Advection in v |  |
| Time-dependent variational principle
---
with RK4
$\mathcal{O}(2D_{W}D^{3}d+D_{W}^{2}D^{2}d^{2})$ | $\mathcal{O}(S^{2}S_{W}^{3}+2S^{4}S^{2}_{W}+S^{4}S_{W})$
| Time-dependent variational principle
---
with exact diagonalization
| $\mathcal{O}(D_{W}^{2}D^{2}d^{2}+D^{4}D_{W}d^{2})$
---
$+(D^{2}d\times D^{2}d)$ eigenvalue problem
| $\mathcal{O}(S^{2}S_{W}^{3}+S^{4}S_{W}^{2}+S^{6}S_{W})$
---
$+(S^{3}\times S^{3})$ eigenvalue problem
Propagation of EM Fields |  |
| Density matrix renormalization group
---
with conjugate gradient descent
$\mathcal{O}(2D_{W}D^{3}d+D_{W}^{2}D^{2}d^{2})$ | N/A
Table 1: Theoretical costs of primary steps in the Vlasov-Maxwell solver for
the tensor train (QTT) geometry and the spine in the comb geometry (QTC-
Spine). In the QTT geometry, all vectors have bond dimension $D$ while
operators have bond dimension $D_{W}$. In the spine of the comb tensor
network, all bonds are of size $S$ and $S_{W}$ for the vectors and operators,
respectively.

##  IV Vlasov-Maxwell Solver

We implement a semi-implicit finite-difference solver of the Vlasov-Maxwell
equations within the quantized tensor network (QTN) format. The distribution
functions for each species and each component of the electromagnetic fields
are represented as QTNs with bond sizes of at most $D$. Likewise, the
operators acting on these fields (e.g., the derivatives) are also represented
in the QTN format (see Ref. Ye and Loureiro (2022) for details).

The distribution functions are updated via explicit time integration of the
Boltzmann equation. We apply second-order operating splitting between the x
and v advection steps. For advection in x, we utilize a split-step semi-
Lagrangian scheme, in which the operator is split further into 1-D advection
steps using a second-order scheme Durran (2010). For advection in v, we
utilize a new time evolution scheme specific to QTNs based on the time-
dependent variational principle (TDVP) Lasser and Su (2022). The
electromagnetic fields are updated using a Crank-Nicolson time integration
scheme, which is solved using an iterative scheme called density matrix
renormalization group (DMRG) Dolgov and Savostyanov (2014).

###  IV.1 Ensuring positivity

One problem with solving the Vlasov equation in the QTN format is that taking
the low-rank approximation of $f$ generally does not preserve important
properties, such as its norm and positivity. The norm can be corrected by
rescaling the distribution function. There are also methods one can employ to
conserve momentum and energy Einkemmer and Jospeh (2021), though these
techniques are not implemented in this work. However, it is difficult to
ensure positivity once in the QTN format.

A simple fix is to define $g$ such that $f=g^{*}g$, and instead evolve $g$
forward in time. Plugging this into Eq. (1), we obtain

$\frac{\partial g_{s}}{\partial
t}+\textbf{v}_{s}\cdot\nabla_{\textbf{x}}g_{s}+\textbf{F}_{s}\cdot\nabla_{\textbf{v}}g_{s}=0$
(9)

and an equivalent equation for $g^{*}$. In the following calculations, we are
actually solving Eq. (9) in order to ensure positivity in our calculations.

###  IV.2 Semi-Lagrangian time evolution

With operator splitting, higher-dimensional advection problems are simplified
to a series of 1-D advection problems of the form

$\frac{\partial\psi}{\partial t}+u\frac{\partial\psi}{\partial x}=0$

where $\psi$ is the field of interest and $u$ is the advection coefficient
that is constant along the $x$ axis. However, $u$ will in general be variable
in all other dimensions of the problem. In the semi-Lagrangian time evolution
scheme, one approximates the field as a set of parcels, computes their exact
trajectories, and then uses the updated parcels to estimate $\psi$ at the next
time step.

Consider a uniform grid, with $x_{j}=j\Delta x$, and let $\phi$ denote a
numerical approximation to field $\psi$; here, $\phi_{j}=\psi(x_{j})$. Since
there are no sources, $\phi$ at the next time step is given by

$\phi(x_{j},t^{n+1})-\phi(\tilde{x}_{j}^{n},t^{n})=0$

where $\tilde{x}_{j}^{n}=x_{j}-u\Delta t$ is the departure point of the
trajectory starting at time $t^{n}$ and arriving at $x_{j}$ at the next time
step $t^{n+1}$. In this work, the value of $\phi(\tilde{x}^{n}_{j})$ is
obtained using a cubic interpolation. For $u>0$, this is given by Durran
(2010)

$\displaystyle\phi(\tilde{x}^{n}_{j},t^{n})=$
$\displaystyle\quad-\frac{\alpha(1-\alpha^{2})}{6}\phi^{n}_{j-p-2}+\frac{\alpha(1+\alpha)(2-\alpha)}{2}\phi^{n}_{j-p-1}$
$\displaystyle\quad+\frac{(1-\alpha^{2})(2-\alpha)}{2}\phi^{n}_{j-p}-\frac{\alpha(1-\alpha)(2-\alpha)}{6}\phi^{n}_{j-p+1}$
(10)

where $p$ is the integer part of $u\Delta t/\Delta x$,
$\alpha=\frac{x_{j-p}-\tilde{x}_{j}^{n}}{\Delta x}$, and superscript $n$
denotes the time step. Thus, the next time step is obtained from
$\phi^{n+1}=A_{SL}\phi^{n}$, where the matrix $A_{SL}$ encodes the above
equation.

Obtaining the QTN-Operator decomposition of $A_{SL}$ can be computationally
expensive. At the moment, we do not have an analytical form for the QTN-
Operator, even for constant $u$, and any expression one does obtain would be
difficult to generalize to higher-dimensional problems with variable $u$. One
can compute the QTN-Operator numerically by performing a series of SVD
operations in an iterative fashion Oseledets and Tyrtyshnikov (2010). However,
the cost of this naive but reliable procedure scales like
$\mathcal{O}(N^{2})$, where $N$ is the total number of grid points along the
dimensions that are not constant. Interpolation algorithms such as TT-Cross
Oseledets and Tyrtyshnikov (2010) often perform well for low-rank tensors,
with computational costs that scale polynomially with rank. However, the
resulting QTN is expected to require moderately large ranks, especially for
larger time steps Kormann (2015).

Fortunately, while semi-Lagrangian time evolution in the QTN format does not
appear to be advantageous for arbitrary problems, it can be relatively
efficient for performing advection of the Vlasov equation along the spatial
dimensions. For example, in the case of advection along the $x$ axis,

$\frac{\partial}{\partial
t}g(\textbf{x},\textbf{v})+v_{x}\frac{\partial}{\partial
x}g(\textbf{x},\textbf{v})=0,$

the operator $A_{SL}$ is effectively 2-dimensional, as it can be written in
the form $A_{SL}=\tilde{A}_{x,v_{x}}\otimes\mathbb{I}_{y,z,v_{y},v_{z}}$ (the
subscripts denote the dimensions in which the operator lives in). It is
observed that the rank of the QTN representation of $\tilde{A}_{x,v_{x}}$
depends on the largest magnitude of integer $p$; so long as moderate time step
and grid resolution is used the bond dimension remains manageable.
Furthermore, if the time step is held constant, the QTN decomposition only
needs to be performed once.

Once in the QTT format, the matrix-vector multiplication is performed using
the density-matrix method tn_ (2023); McCulloch (2007). Details of this
algorithm for the comb geometry are given in Supplementary Information (SI),
and the computational costs are listed in Table 1.

###  IV.3 Time-dependent variational principle

Thus far, prior works with QTNs utilize classical explicit time integrators
(such as fourth-order Runge-Kutta (RK4), MacCormack). Similar to the semi-
Lagrangian scheme discussed above, most explicit schemes only require matrix-
vector multiplications and additions, and thus can easily be translated into
the QTN format. These time integration schemes are also global schemes, since
they involve manipulating the entire QTN state at once.

We can instead consider a local time evolution scheme, in which individual
tensors in the QTN are propagated forward in time in a sequential fashion. In
quantum simulation, one such algorithm is the time-dependent variational
principle (TDVP) scheme Haegeman  _et al._ (2016). This scheme is also closely
related to the dynamical low-rank algorithm Khoromskij  _et al._ (2012);
Einkemmer and Lubich (2018), which has been used for solving PDEs in the not
quantized tensor train format (only utilizing low-rank tensor decompositions
between but not within each problem dimension). More recently, the same
underlying principle has been used in a quantum algorithm for solving the
Schrodinger-Poisson equation Cappelli  _et al._ (2023).

There are very fundamental differences between global and local time
evolution. In global time evolution, the state is first propagated in time and
then projected onto the manifold for QTT of rank $D$. In contrast, in local
schemes, the dynamics themselves are approximated by projection onto that
manifold.

TDVP utilizes the Dirac-Frenkel variational principle, which states that the
dynamics of a unitary system at any given point in time lies in the tangent
space of the current state Frenkel (1934); Lasser and Su (2022). By treating
each tensor as a variational parameter, one can obtain the effective dynamics
for each tensor (given by anti-Hermitian operator $A_{\text{eff}}^{(i)}$ for
the $i^{\text{th}}$ tensor) and sequentially evolve them forward in time. More
specifically, if the unitary dynamics acting on state $\psi$ is described by
anti-Hermitian operator $A$, the equation of motion for the $i^{\text{th}}$
tensor $M_{i}$ is

$\displaystyle\frac{\partial}{\partial t}M_{i}$
$\displaystyle=A_{\text{eff}}M_{i}\,,$ (11)

with

$\displaystyle
A_{\text{eff}}=\partial_{M_{i}}\,\partial_{M_{i}^{*}}\langle\psi(M_{1},...,M_{L})\,|\,A\,|\,\psi(M_{1},...M_{L})\rangle$
(12)

assuming that the QTT is canonicalized such that the orthogonality center is
located at tensor $i$. As a result, the original dynamics is decomposed into
roughly $2KL$ equations of motion of size $~{}\mathcal{O}(D^{2})$. (The factor
of two arises from the algorithm when sweeping across all tensors.) Eq. (11)
is solved using RK4, so the computational cost scales like
$\mathcal{O}(D^{3}D_{W}d+D^{2}D_{W}^{2}d^{2})$, where $D_{W}$ is the bond
dimension of the QTT-O representation for time evolution operator $A$. We
refer the reader to the SI and Ref. Paeckel  _et al._ , 2019 for algorithmic
details.

The time evolution scheme presented above is the single-site variant of TDVP.
By construction, the time-evolved state will always have the same bond
dimension as the initial state. As a result, the algorithm is cheaper than
global methods since it avoids needing to compress the state from an enlarged
bond dimension. However, this feature can be problematic when the time evolved
state requires larger bond dimension, particularly when the initial bond
dimension is very small. To counter this, one would typically use the slightly
more expensive two-site variant in which two neighboring tensors are updated
simultaneously so the QTN bond dimension can be adjusted as needed. However,
in our solver, single-site TDVP appears to be sufficient because it is used in
conjunction with global time evolution schemes that already increase QTN rank.
For both variants, we found it imperative to use a global time evolution
scheme such as RK4 for the first time step.

Particularly interesting is that even when Eq. (11) is solved approximately
using RK4, we observe that the local methods allow one to use larger time
steps than what global methods are limited to (e.g., by the Courant-Friedrich-
Levy (CFL) constraint Courant  _et al._ (1967)). In the following example, it
appears that the time steps of the local algorithm can be larger by at least a
factor of four. This is perhaps possible because global information is
incorporated into the effective operator for each individual tensor. The
ability to bypass CFL is crucial for QTT finite difference solvers, as it
allows one to use higher spatial resolution (for improved accuracy) without
requiring one to significantly increase temporal resolution.

####  IV.3.1 Advection with time varying electric field

We demonstrate TDVP time evolution with the advection of a 0D2V uniform
electron test distribution in the presence of constant background magnetic
field $B_{0,z}\hat{z}$, and uniform but time-varying electric field
$\textbf{E}_{x}(t)=E_{0}\cos(\omega t)\hat{x}$. The analytical solution to the
dynamics results in the drift velocities

$\displaystyle v_{x}(t)=-\omega_{x}(t)-v_{x,0}\cos(t)-v_{y,0}\sin(t)$ (13)
$\displaystyle v_{y}(t)=-\omega_{y}(t)+v_{x,0}\sin(t)-v_{y,0}\cos(t)$ (14)

with

$\displaystyle\omega_{x}(t)$
$\displaystyle=\frac{E_{0}}{1-\omega^{2}}\left(\sin(t)-\omega\sin(\omega
t)\right)$ (15) $\displaystyle\omega_{y}(t)$
$\displaystyle=\frac{E_{0}}{1-\omega^{2}}\left(\cos(t)-\cos(\omega t)\right)$
(16)

for $\omega\neq 1$. The drift velocity in the $z$ direction remains constant.

We perform time evolution using a second-order semi-Lagrangian scheme, fourth-
order Runge-Kutta (RK4), and TDVP. The calculations use $2^{L}$ grid points
along the $v_{x}$ and $v_{y}$ axes (with bounds $[-12v_{th},12v_{th})$). The
test electron distribution is initialized with drift velocities
$v_{x,0}=v_{y,0}=0$ and thermal velocity of 1, and is evolved until a time of
100 using time step $\Delta t$. The calculations were performed using a
sequential layout with data along both axes encoded via a binary mapping. Both
second-order and fourth-order centered finite difference stencils are
considered.

Figure 2: Rms errors in the dynamics for advection with time-varying electric
field. Test parameters are set to be $\omega=0.4567$, $E_{0}=0.9$, and
$B_{0,z}=1$. Calculations are performed with bond dimension $D$, grid
resolution $2^{L}$ along each dimension, and time step $\Delta t$. (top)
Results for varied $D$ and fixed $L=8$ and $\Delta t=0.025$. (middle) Results
for varied $L$ and fixed $D=16$ and $\Delta t=0.1$. (bottom) Results for
varied $\Delta t$ and fixed $L=8$ and $D=16$. The different plots along each
row are labeled with the time-evolution scheme used: semi-Lagrangian (SL),
fourth-order Runge-Kutta (RK4), and single-site TDVP. Errors with respect to
the variance (i.e., thermal velocity squared, green), drift velocity (blue),
and smoothness of the Maxwellian distribution (Eq. (17); red) are plotted in
different colors. The two symbols denote calculations performed with a second-
order or fourth-order centered finite-difference stencil. By construction, the
stencil plays no role in SL time evolution. Missing data is a result of
numerical instabilities.

At each time step, the distribution function $|g(v_{x},v_{y})|^{2}$ obtained
after time evolution is fit to a Maxwellian $\mathcal{M}$ and the errors in
the variance ($\epsilon_{\sigma}$) and drift velocity ($\epsilon_{v}$) are
computed as the norm of the differences of the fit with respect to the
theoretical values along the two dimensions. The “error in the smoothness” is
quantified as

$\epsilon_{\nabla}=\frac{1}{2^{L}}\left\lVert\left(\frac{d}{dv_{x}}+\frac{d}{dv_{y}}\right)\left(|g(v_{x},v_{y};t)|^{2}-\mathcal{M}_{\text{fit}}\right)\right\rVert$
(17)

where the derivatives are computed with a first-order finite difference
stencil. The total error of the simulation is the root-mean-square error,
$\lVert\epsilon(t)\rVert/\sqrt{N_{t}}$, where $N_{t}$ is the total number of
time steps.

The results are shown in Fig. 2. Plots in the top row show that a bond
dimension of $D=16$ is sufficient for these calculations, and the error is not
reduced with increased bond dimension. In the single-site TDVP calculations,
the bond dimension throughout the simulation is fixed at the bond dimension of
the initial state, which is about $D=10$. Plots in the middle row show how
error scales with increasing grid resolution, for fixed bond dimension $D=16$
and time step $\Delta t=0.1$. This time step is too large for RK4 with
resolutions greater than $L=6$, leading to numerical instabilities. In the SL
case, increasing resolution most notably improves the error in the variance
(temperature) of the Maxwellian, as well of the smoothness of the resulting
distribution function. For TDVP, increasing resolution also reduces errors,
though at higher grid resolutions another source of error seems to appear. The
bottom row compares errors for calculations with $D=16$, $L=8$ for different
time steps. With SL time evolution, error in the drift velocity decreases with
reduced time step, but the error in variance and smoothness increases. This
suggests that the error arises from the projection of the Lagrangian dynamics
onto the grid at each time step. With TDVP, for the calculation with the
fourth-order finite-difference stencil, performance improves with reduced time
step, as is expected since the sequential iterative updates have an associated
commutation error. However, with the second-order stencil, the error appears
to be dominated by lack of resolution since reducing time step does not reduce
error.

In summary, these results show that both SL time evolution and TDVP allow one
to take larger time steps than with conventional grid-based methods like RK4.
However, note that when TDVP is performed with RK4 as the local time
integrator (as is done here), there still exists a maximum time-step above
which time evolution becomes unstable. In general, the error in the thermal
velocity is smaller with TDVP than SL time evolution (with which the
Maxwellian tends to broaden). However, the distribution function from SL tends
to be smoother than the results from TDVP. This may become problematic for
long-time time integration using TDVP.

###  IV.4 Implicit solver for Maxwell’s equations

The electric and magnetic fields are propagated forward in time implicitly
using Crank-Nicolson. Given the fields at time step $n$, the updated fields at
the next time step are obtained by solving

$\displaystyle\begin{bmatrix}\frac{1}{\Delta
t}\mathbb{I}&-\frac{c^{2}}{2}\nabla\times\\\\[5.69046pt]
\nabla\times&\frac{1}{\Delta
t}\mathbb{I}\end{bmatrix}\begin{bmatrix}\textbf{E}^{n+1}\\\\[5.69046pt]
\textbf{B}^{n+1}\end{bmatrix}=$ $\displaystyle\hskip
14.22636pt\begin{bmatrix}\frac{1}{\Delta
t}\mathbb{I}&\frac{c^{2}}{2}\nabla\times\\\\[5.69046pt]
-\nabla\times&\frac{1}{\Delta
t}\mathbb{I}\end{bmatrix}\begin{bmatrix}\textbf{E}^{n}\\\\[5.69046pt]
\textbf{B}^{n}\end{bmatrix}+\begin{bmatrix}c^{2}\nabla\times\textbf{B}_{0}-\frac{1}{\varepsilon_{0}}\textbf{J}^{n}\\\\[5.69046pt]
-\nabla\times\textbf{E}_{0}\end{bmatrix}$ (18)

where $\textbf{E}_{0}$ and $\textbf{B}_{0}$ are constant background
contributions to the electric and magnetic fields, and $\textbf{J}^{n}$ is the
current density measured from the distribution functions at the
$n^{\text{th}}$ time step.

Linear equations of the form $Ax=b$ can be solved approximately in the QTT
format using the density matrix renormalization group (DMRG) algorithm
Oseledets and Dolgov (2012); Dolgov and Savostyanov (2014); Schollwöck (2011);
Gourianov  _et al._ (2022), in which each tensor (or pairs of neighboring
tensors) are optimized sequentially to minimize the cost function
$x^{\dagger}Ax-2x^{\dagger}b$. (Details are provided in the SI.) If the local
optimizations are performed using conjugate gradient descent (or conjugate
gradients squared Sonneveld (1989) for non-Hermitian $A$), then the algorithm
primarily consists of tensor contractions and the cost of the algorithm scales
like $\mathcal{O}((D_{x}^{2}D_{b}+D_{b}D_{x}^{2})D_{A}^{2})$, where $D_{x}$
and $D_{b}$ are the bond dimensions of QTT-States $x$ and $b$, and $D_{A}$ is
the bond dimension of QTT-Operator $A$. For the operator in Eq. (18), $D_{A}$
is moderate and independent of system size, since $A$ can be built directly
from differential operators and the calculations are performed on a uniform
Cartesian grid.

If the QTT bond dimension is sufficiently large and the time-step is small
enough such that matrix $A$ is well-conditioned, then DMRG converges quickly
to the expected solution. In practice, DMRG often converges to some value
larger than the tolerated error threshold (of $10^{-6}$). In this case, we
split a single time step into at most 4 smaller equally sized time steps; if
DMRG shows better convergence at each step, then we continue to solve
Maxwell’s equations with these reduced time steps. Otherwise, we simply accept
that the QTT bond dimension is likely not large enough and continue with the
simulation as is. Certainly, a more advanced implicit time evolution scheme
should be considered in the future.

Note that the DMRG algorithm can be extended to more general QTN geometries.
However, in this work we are only solving Maxwell’s equation in two
dimensions, and the QTC with two branches is equivalent to a QTT. As such, we
do not need to modify the original algorithm.

##  V Results

In the following section, we evaluate the performance of the proposed QTN
methods on two common test problems. Our code is built off the tensor network
software package quimb Gray (2018).

###  V.1 Orszag-Tang

(a) $D=32$, $t=21\Omega_{c,p}^{-1}$ |  (b) $D=64$, $t=21\Omega_{c,p}^{-1}$ |  (c) $D_{f}=32,\,D_{EM}=128$, $t=9\Omega_{c,p}^{-1}$
---|---|---

(a)

Figure 3: Plots of the out-of-plane current for the Orszag-Tang vortex at the
specified time for (a) $D=32$, (b) $D=64$, and (c) bond dimension $D_{f}=32$
for the distribution functions and bond dimension $D_{EM}=128$ for the
electromagnetic fields. Calculations were performed with a time step of
$0.001\Omega_{c,p}^{-1}.$

We consider a (fully nonlinear and kinetic) 2D3V calculation of the Orszag-
Tang vortex Orszag and Tang (1979) with protons and electrons. The initial
conditions are Juno  _et al._ (2018)

B
$\displaystyle=-B_{0}\sin\left(k_{y}y\right)\hat{x}+B_{0}\sin\left(2k_{x}x\right)\hat{y}+B_{g}\hat{z}$
(19) $\displaystyle\textbf{u}_{p}$
$\displaystyle=-u_{0,p}\sin\left(k_{y}y\right)\hat{x}+u_{0,p}\sin\left(k_{x}x\right)\hat{y}$
(20) J
$\displaystyle=en_{0}(\textbf{u}_{p}-\textbf{u}_{e})=\frac{1}{\mu_{0}}\nabla\times\textbf{B}$
(21) E $\displaystyle=-\textbf{u}_{e}\times\textbf{B}$ (22)

where $k_{x}=k_{y}=2\pi/L_{\text{box}}$, with $L_{\text{box}}$ being the
simulation domain size in real space. $B_{g}$ is the out-of-plane guide field.
The amplitudes of the perturbations are set to $B_{0}=0.2B_{g}$ and
$u_{0,p}=0.2v_{A,p}$, where $v_{A,p}=B_{g}/\sqrt{\mu_{0}n_{0}m_{p}}$ is the
Alfvén velocity of the protons. The initial velocities ($\textbf{u}_{s}$) are
included as initial drift velocities of the Maxwellian distribution,

$\displaystyle f_{s}(\textbf{v})=\left(2\pi
v_{th,s}^{2}\right)^{-3/2}\exp\left(-\frac{(\textbf{v}-\textbf{u}_{s})^{2}}{2v_{th,s}^{2}}\right),$
(23)

where $v_{th,s}^{2}=k_{B}T_{s}/m_{s}$ is the thermal velocity of species $s$,
with $k_{B}$ the Boltzmann constant, $m_{s}$ the mass, and $T_{s}$ the
temperature.

The plasma is initialized with density $n_{0}=1$, mass ratio $m_{p}/m_{e}=25$,
temperature ratio of $T_{p}/T_{e}=1$, proton plasma beta
$\beta_{p}=2v^{2}_{th,p}/v^{2}_{A,p}=0.08$ and speed of light $c=50v_{A,p}$.

The simulation is performed with $L_{\text{box}}=20.48d_{p}$, where
$d_{p}=\sqrt{m_{p}/\mu_{0}n_{0}e^{2}}$ is the proton skin depth, with $e$
being the elemental charge. The limits of the grid in velocity space for
species $s$ are $\pm 7v_{th,s}$. Periodic boundary conditions are used in all
dimensions, including the velocity dimensions, so that the time evolution
operator is unitary and TDVP can be used. The calculations are performed with
$2^{9}$ grid points along each spatial axis, and $2^{6}$ grid points along
each axis in velocity space. With this grid resolution, the grid-spacing in
real space is 0.04$d_{p}$, so the electron skin depth is well-resolved. We
utilize the QTC ansatz with tensors corresponding to each dimension on its own
branch, ordered as $(x,y,v_{x},v_{y},v_{z})$, and the data is encoded with
binary mapping such that fine grid tensors are attached to the spine for all
dimensions. At each time step, the proton and electron distribution functions,
as well as each component of the electric and magnetic fields, are compressed
to bond dimension $D$.

Fig. 3 shows the out-of-plane current density at time $21\Omega_{c,p}^{-1}$
obtained from simulations with bond dimension (a) $D=32$ and (b) $D=64$.
Qualitatively, both calculations exhibit the expected dynamics. However, the
features in the $D=32$ calculation appear less sharp than in the $D=64$ case,
and there is significantly more numerical noise. Unfortunately, the numerical
noise will continue to grow for both calculations. Visual inspection suggests
that the noise arises from the Maxwell equation solver in the electric field
(see the SI for these plots). However, when increasing the bond dimension of
the QTTs representing the electric and magnetic fields to $128$, the noise
actually worsens and one is limited to even shorter simulation times (Fig.
3(c)). This suggests that the low-rank approximation provides some noise
mitigation; with less compression, the noise grows unchecked.

Random and unstructured functions are not compressible in the QTT format. As
such, numerical noise will likely reduce the compressibility of the original
function. Unfortunately, it is unlikely one will obtain the original function
when making low-rank approximation on a noisy function. Instead, it is more
likely that we obtain low-rank noise that exhibits some structure. Given the
observation that the features in the $D=32$ calculation are broader than in
the $D=64$ case, it appears that after the noise accumulates over time, the
QTT is only able to capture features at an “effective grid resolution”
determined by $D$, i.e., the grid resolution for which the specified $D$ would
yield no compression. For example, in these calculations, $D=32$ for the
electromagnetic fields (which are 2D) would correspond to a grid resolution of
$0.64d_{p}$, while $D=64$ would correspond to a grid resolution of
$0.32d_{p}$. These length scales also appear as special scales in the spectra
(see Section V.A of the SI).

###  V.2 GEM Reconnection

(a) $D=32$

(b) $D=64$

Figure 4: Out-of-plane current density ($J_{z}$) for calculations with bond
dimension (a) $D=32$ and (b) $D=64$. Calculations were performed with time
step $\Delta t=0.005\Omega_{c,p}^{-1}$. Time is in units of
$\Omega_{c,p}^{-1}.$

The next test we consider is the classic GEM reconnection problem Birn and
Hesse (2001) in 2D3V, again reformulated for fully kinetic simulation Schmitz
and Grauer (2006); Allmann-Rahn  _et al._ (2022). Some additional
modifications are made such that periodic boundary conditions can be used.

In the traditional set-up, the initial configuration for each species $s$ is

$\displaystyle f_{s}=f_{0,s}+f_{b,s}$ (24)

where

$\displaystyle f_{b,s}=n_{b}(2\pi
v_{th,s}^{2})^{-3/2}\exp(-\textbf{v}^{2}/2v_{th,s}^{2})$ (25)

is a uniform background distribution and $f_{0,s}$ is the Harris sheet
equilibrium Harris (1962). Given the magnetic field

$\displaystyle\textbf{B}(y)=B_{0}\tanh\left(\frac{y}{\lambda}\right)\hat{x},$
(26)

where $\lambda$ is the thickness of the sheet, the proton/electron
distribution is

$\displaystyle f_{0,s}=\frac{n_{0}(y)}{(2\pi
v_{th,s}^{2})^{3/2}}\exp\left(-\frac{v_{x}^{2}+v_{y}^{2}+(v_{z}-v_{0,s})^{2}}{2{v^{2}_{th,s}}}\right),$
(27)

where

$\displaystyle
n_{0}(y)=n_{0}\operatorname{sech}^{2}\left(\frac{y}{\lambda}\right).$ (28)

The drift velocities are constrained by the equilibrium conditions,

$\displaystyle\frac{B_{0}^{2}}{2\mu_{0}}=n_{0}(T_{p}+T_{e}),$ (29)
$\displaystyle\lambda=\frac{2}{eB_{0}}\frac{T_{p}+T_{e}}{v_{0,p}-v_{0,e}},$
(30) $\displaystyle\frac{v_{0,e}}{v_{0,p}}=-\frac{T_{e}}{T_{p}}.$ (31)

The initial perturbation in the magnetic field is

$\displaystyle B_{x}$
$\displaystyle=-\frac{A_{0}k_{y}}{2}\cos(k_{x}x)\,\sin(k_{y}y/2)\,,$ (32)
$\displaystyle B_{y}$
$\displaystyle=A_{0}k_{x}\sin(k_{x}x)\,\cos(k_{y}y/2)\,.$ (33)

To allow for periodic boundary conditions, the simulation domain (suppose
previously $\left[-L_{\text{box},x}/2,L_{\text{box},x}/2\right)$ and
$\left[-L_{\text{box},y}/2,L_{\text{box},y}/2\right)$) is doubled in $y$, and
an inverted Harris sheet placed at the simulation boundary:

$\displaystyle B_{x}(y)=$ $\displaystyle
B_{0}\tanh\left(\frac{y}{\lambda}\right)$ $\displaystyle-
B_{0}\tanh\left(\frac{y+L_{\text{box},y}}{\lambda}\right)$ $\displaystyle-
B_{0}\tanh\left(\frac{y-L_{\text{box},y}}{\lambda}\right)\,,$ (34)

and

$\displaystyle n(y)=$ $\displaystyle
n_{0}\operatorname{sech}((y+L_{\text{box},y})/\lambda^{2})$
$\displaystyle+n_{0}\operatorname{sech}((y-L_{\text{box},y})/\lambda^{2})\,.$
(35)

The perturbation in the magnetic field remains the same.

Simulations are performed with $L_{\text{box},x}=25.6d_{p}$,
$L_{\text{box},y}=12.8d_{p}$, and a Harris sheet thickness of
$\lambda=0.5d_{p}$. The plasma is defined to have a mass ratio of
$m_{p}/m_{e}=25$, temperature ratio of $T_{p}/T_{e}=5$, proton plasma beta
$\beta_{p}=5/6$, and electron Alfvén velocity $v_{A,e}=0.25c$. The strength of
the perturbation is $A_{0}=0.1B_{0}d_{p}$. The background density is set to
$n_{b}=0.2n_{0}$. The limits of the grid in velocity space are $\pm
7.75v_{th,p}$ for the protons. For the electrons, the limits are $\pm
7v_{th,e}$ in the $v_{x}$ and $v_{y}$ directions, and $\pm 14v_{th,e}$ in the
$v_{z}$ direction. The resolution of the calculations are $2^{8}\times
2^{9}\times 2^{7}\times 2^{7}\times 2^{8}$ along the $x,\,y,\ v_{x},\,v_{y},$
and $v_{z}$ axes respectively. We again utilize the QTC ansatz with tensors
corresponding to each dimension on its own branch. In contrast to above, the
data is encoded with the mirror (two’s complement) mapping for all dimensions
because it showed marginally better results. Additionally, the bond dimension
for the electric and magnetic fields are set to 128 so that the DMRG solver
for Maxwell’s equations might converge more quickly.

(a)  Figure 5: Reconnected magnetic flux computed with time step $\Delta
t=0.005\Omega_{c,p}^{-1}$ and bond dimension $D=32$ (blue) and $D=64$
(orange). Time is in units of $\Omega_{c,p}^{-1}$.

Fig. 4 shows the current density and Fig. 5 plots the reconnected magnetic
flux $\Psi=\int_{0}^{L_{\text{box},x}/2}B_{y}(x,y=0)\,dx$ for calculations
with bond dimensions $D=32$ and $D=64$. Additional results are presented in
Section V.B of the SI. Again, for both bond dimensions considered here, we
roughly obtain the expected reconnected flux over time Schmitz and Grauer
(2006). However, the dynamics appear sensitive to numerical noise, as
evidenced by unexpected secondary island formation in the Harris sheet for the
$D=32$ case.

##  VI Discussion

The above results are a promising first demonstration that high-dimensional
simulations of the Vlasov-Maxwell equations can be performed efficiently with
the QTC ansatz. Our calculations use $N=2^{36}$ total grid points for the
Orszag-Tang vortex, and $N=2^{39}$ total grid points for the GEM reconnection
problem. However, we are able to observe the expected phenomenological results
satisfactorily with a maximum bond dimension of just $D=64$. Despite the
computational overhead in the QTN algorithm, ignoring constant factors, the
$\mathcal{O}(D^{4})$ cost per time step is a fraction of the $\mathcal{O}(N)$
scaling for traditional methods, and the simulations can be run on a single
compute node. (Tighter bounds can be obtained by distinguishing between the
virtual bonds along the spine and along the branches). Furthermore, by using
the local time evolution scheme TDVP, one can use a larger time step (observed
to be larger by about a factor of 4) than allowed in global time evolution
schemes. In addition to computational cost, memory storage costs are also
reduced, from $\mathcal{O}(N)$ to $\mathcal{O}(D^{3})$.

The von Neumann or bipartite entanglement entropy (EE) Schollwöck (2011) is
one measure of efficiency of the QTC representation. If considering the
$j^{\text{th}}$ bond in the tensor network, it is given by

$EE^{(j)}=-\sum_{i}\tilde{\sigma}_{i}^{2}\log_{2}(\tilde{\sigma}_{i}^{2})$
(36)

where $\tilde{\sigma}$ are the singular values at bond $j$ when the tensor
network’s center of center of orthogonality is located at bond $j$. (The
singular values are normalized such that $\sum_{i}\tilde{\sigma}_{i}^{2}=1$.)
We compute the EE at each bond in the QTC for the ion distribution, electron
distribution, electric fields, and magnetic fields. Results for the Orszag-
Tang vortex are shown in Fig. 6, and results for the reconnection problem are
shown in Fig. 7. Interestingly, in both cases, the EE of the ion and electron
distributions remain relatively small, and the EE of the electric and magnetic
fields tend to grow. Some of that is likely due to noise, as suggested by the
sudden increase in EE and the visible increase of noise in the fields
themselves (see the SI for these figures). Most of the noise appears to be
injected through the electric field. This suggests that the issue is likely
not due to insufficient bond dimension but rather numerical noise injected
through various stages of the solver. This may be because our current QTT
format solvers (i.e., TDVP and DMRG) are prone to introducing noise, but it
also may be due to the simulation parameters used—the time step may have been
too large for the Crank-Nicolson scheme, and the grid resolution used still is
not fine enough to resolve the electron Debye length. Fortunately, it should
be possible to remedy both types of issues in the future.

(a)

(b)

Figure 6: Entanglement entropy at each bond for the Orszag-Tang problem. Bond
number is ordered by bonds in the spine, and then bonds in each branch.

(a)

(b)

Figure 7: Entanglement entropy at each bond for the reconnection problem. Bond
number is ordered by bonds in the spine, and then bonds in each branch.

As with most local optimization schemes, the DMRG solver used for the implicit
time evolution of Maxwell’s equations is prone to slow convergence and only
finding local minima, particularly when the bond dimension is small. Our
current implementation is very simple; using a preconditioner and a better
initial guess will likely improve performance. One might also benefit from
adding additional constraints, such as one to preferentially find smooth
results, to help improve convergence. Whether or not these constraints can be
also be used to improve numerical stability of the implicit solver as well is
also worth more in-depth investigation. Alternatively, one could include a
filter or artificial dissipation that removes grid-scale noise.

In addition to improving our solver for Maxwell’s equations, we also need to
improve our understanding of the TDVP time evolution scheme. In particular,
while single-site TDVP time evolution is used in this work with promising
results, local tangent-space time evolution schemes need be examined further,
since they have their limitations. Most notably, TDVP (both 1-site and 2-site)
is known to sometimes yield incorrect dynamics Yang and White (2020), even as
one reduces the time step. Understanding the lack of convergence with respect
to time step and grid resolution must be addressed in order to trust results
for problems we do not know the solution to. Additionally, when using TDVP
with RK4, there still is a time-step constraint, but we currently do not have
an estimate for it. We also need to extend TDVP to non-Hermitian systems to
account for non-periodic boundary conditions, so that a greater variety of
problems can be considered. It may be possible to improve our basic
implementation to better ensure conservation laws. For example, for simulation
of closed quantum systems, the TDVP algorithm naturally conserves system
energy and system norm, because the Hamiltonian describes both system energy
and dynamics. Though this is not true in the case of the Vlasov equation,
there may be an alternative formulation that allows one to take advantage of
these conservation properties. There also may be ways to impose additional
conservation properties, as is done in the dynamical low-rank literature
Einkemmer and Lubich (2019); Einkemmer  _et al._ (2020); Einkemmer and Jospeh
(2021).

As of now, our code runs serially on a single CPU node, with the exception of
the internal threading performed by Python’s numpy package. As such, while the
amount of resources required to perform and store these calculations is much
less than conventional methods, the overall run-time is likely longer than
that of an optimized solver. While our code could likely be optimized, a more
fundamental problem is that the algorithm for QTTs is inherently serial, thus
making parallel calculations unattractive. The QTC geometry may be more
amenable for parallelization, with algorithms modified to allow one to work
with each branch independently. However, computations on the spine generally
have worse scaling than the branches, and would face the same challenges for
parallelization as QTTs. The parallelized algorithms likely will also be less
efficient than serial algorithms, requiring repetitive calculations. For
example, instead of computing the necessary information for just one branch
and then building the information for the remaining branches as one proceeds
through the algorithm, one may need to compute the required information for
each branch all at once. However, given the low computational costs of the
serial algorithm for modest bond dimension, increasing the cost linearly with
problem dimension for parallelization is affordable.

One might argue that the observed efficiency and compressibility of the QTN
representation is due to the inefficiency of the grid-based representation,
and one could instead use a spectral or finite element basis. Ref. Juno  _et
al._ , 2018 reports a discontinuous Galerkin finite element calculation of the
Orszag-Tang vortex with similar parameters. They use a $40^{2}\times 8^{3}$
grid in 2D3V and a second-order polynomial basis, resulting in $N_{p}=112$
degrees of freedom per grid cell. Thus, the distribution function is
represented using $NN_{p}$ ($9\times 10^{7}$) degrees of freedom (d.o.f.),
which is greater than the approximate upper bound of $5D^{3}$ for the QTC
ansatz ($1.3\times 10^{6}$ d.o.f. for $D=64$). By this rough comparison, the
QTC ansatz indeed provides a smaller but approximate representation of the
distribution function. For reference, newer finite element simulations on a
$112^{2}\times 16^{3}$ grid up to a time of $75\Omega_{c,p}^{-1}$ used 10k
node-hours Juno  _et al._ (2018); Hakim  _et al._ (2020); Hakim and Juno
(2020).

In comparison to non-quantized low-rank methods, the proposed QTN method
exhibits a much more modest reduction in computational resources for the
problems discussed above. One advantage of the QTN methods is the ability to
use time-steps larger than the CFL limit because of the local time evolution
scheme, without having to perform semi-Lagrangian time evolution. Furthermore,
QTNs are more likely to demonstrate a bigger advantage for problems requiring
orders of magnitude more grid points along each dimension, as will be
necessary for more realistic problems (e.g., with the proper mass ratio). A
more thorough investigation on the efficiency of the QTN format for different
simulation parameters needs to be performed, and testing the proposed
algorithm for more complex problems, potentially in different coordinate
systems, is certainly worth pursuing. The QTN format is also not limited to a
real-space, grid-based interpretation. Since the convolution operation can be
performed efficiently within the QTN format Kazeev  _et al._ (2013b), one can
also implement an efficient QTN solver with spectral bases. Existing work also
discusses finite-element solvers in the QTN format Kazeev  _et al._ (2017);
Kazeev and Schwab (2018); Kazeev  _et al._ (2022); Kornev  _et al._ (2023).
Performing a careful comparison to understand the benefits and caveats of each
representation would be of practical interest and may improve our
understanding of QTNs more generally.

There are many immediate and practical uses for QTN-based solvers. Because the
low-rank approximations in the QTN format are solely numerical, one is able to
solve the Vlasov equation at reduced cost without relying on assumptions that
would limit the range of problems that can be studied. As such, one could
perform a rough parameter sweep using the proposed solver before identifying a
specific regime in which more accurate calculations are needed. The QTN
solvers can also help with generating approximate but physically accurate
training data for machine learning methods.

Lastly, we also hope that the QTN representation can act as another tool for
understanding turbulence. Drawing from ideas in quantum information theory,
one can better understand information flow between different length scales.
Additionally, by connecting QTNs to other concepts like wavelets Oseledets and
Tyrtyshnikov (2011) and renormalization group, one may be able to build a
consistent picture between different interpretations.

##  VII Conclusion

Solving the Vlasov equation using traditional algorithms is extremely resource
intensive, so an alternative approach that offers reduced computational cost
but still captures kinetic physics is very much desired. In this work, we
implemented a grid-based Vlasov-Maxwell solver utilizing the quantized tensor
network (QTN) format. We utilize a semi-Lagrangian time integrator, as well as
a QTN time evolution scheme that allows one to overcome the CFL time-step
limit that plagues high resolution calculations. The cost of the algorithm
scales roughly like $\mathcal{O}(D^{4})$, where $D$ is the rank or bond
dimension of the QTN, and determines the accuracy of the calculation. In the
test cases considered, we find that QTNs of $D=64$ are sufficient to capture
the expected dynamics with reasonable accuracy, even though the simulation is
performed on grids of around $2^{39}$ grid points in total. Furthermore, the
QTN algorithm allows one to take larger time steps than traditional methods,
further reducing the overall cost. These results suggest that QTN-based
algorithms may enable higher resolution simulations of the Vlasov equation
than otherwise feasible.

##  VIII Code Availability

Code can be made available upon request.

###### Acknowledgements.

EY thanks Jimmy Juno for discussions regarding the cost of the discrete-
Galerkin finite element Vlasov solver in Gkeyll. The authors were supported by
awards DE-SC0020264 and DE-SC0022012 from the Department of Energy.

## References

  * Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
  * Vidal (2003) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
  * White (1993) S. R. White, Phys. Rev. B 48, 10345 (1993).
  * Chan and Head-Gordon (2002) G. K.-L. Chan and M. Head-Gordon, J. Chem. Phys. 116, 4462 (2002).
  * White (2005) S. R. White, Phys. Rev. B 72, 180403(R) (2005).
  * Khoromskij (2011) B. N. Khoromskij, Constr. Approx. 34, 257 (2011).
  * Dolgov  _et al._ (2012) S. V. Dolgov, B. N. Khoromskij, and I. V. Oseledets, SIAM J. Sci. Comput. 34, A3016 (2012).
  * Oseledets (2009) I. V. Oseledets, Dokl. Mathematics 80, 653 (2009).
  * Oseledets (2010) I. V. Oseledets, SIAM J. Sci. Comput. 31, 2130 (2010).
  * Kazeev and Khoromskij (2012) V. A. Kazeev and B. N. Khoromskij, SIAM J. Matrix Anal. A. 33, 742 (2012).
  * Kazeev  _et al._ (2013a) V. Kazeev, O. Reichmann, and C. Schwab, Linear Algebra. Appl. 438, 4024 (2013a).
  * Gourianov  _et al._ (2022) N. Gourianov, M. Lubasch, S. Dolgov, _et al._ , Nature Comput. Sci. 2, 30 (2022).
  * Kiffner and Jaksch (2023) M. Kiffner and D. Jaksch, “Tensor network reduced order models for wall-bounded flows,”  (2023), arXiv:2303.03010 [physics, physics:quant-ph], arXiv:2303.03010 [physics.flu-dyn] .
  * Kornev  _et al._ (2023) E. Kornev, S. Dolgov, K. Pinto, M. Pflitsch, M. Perelshtein, and A. Melnikov, “Numerical solution of the incompressible Navier-Stokes equations for chemical mixers via quantum-inspired Tensor Train Finite Element Method,”  (2023), arXiv:2305.10784 [physics.flu-dyn] .
  * Peddinti  _et al._ (2023) R. D. Peddinti, S. Pisoni, A. Marini, P. Lott, H. Argentieri, E. Tiunov, and L. Aolita, “Complete quantum-inspired framework for computational fluid dynamics,”  (2023), arXiv:2308.12972 [physics.flu-dyn] .
  * Kormann (2015) K. Kormann, SIAM J. Sci. Comput. 37, B613 (2015).
  * Einkemmer and Lubich (2018) L. Einkemmer and C. Lubich, SIAM J. Sci. Comput. 40, B1330 (2018).
  * Allmann-Rahn  _et al._ (2022) F. Allmann-Rahn, R. Grauer, and K. Kormann, J. Comput. Phys. 469, 111562 (2022).
  * Einkemmer (2023) L. Einkemmer, “Accelerating the simulation of kinetic shear Alfvén waves with a dynamical low-rank approximation,”  (2023), arXiv:2306.17526 [physics.comp-ph] .
  * Ye and Loureiro (2022) E. Ye and N. F. G. Loureiro, Phys. Rev. E 106, 035208 (2022).
  * Dolgov (2014) S. Dolgov,  _Tensor product methods in numerical simulation of high-dimensional dynamical problems_ , Ph.D. thesis, Leipzig University (2014).
  * Ripoll (2021) J. J. G. Ripoll, Quantum 5, 431 (2021).
  * Lubasch  _et al._ (2020) M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, Phys. Rev. A 101, 010301 (2020).
  * Jaksch  _et al._ (2023) D. Jaksch, P. Givi, A. J. Daley, and T. Rung, AIAA J. 61, 1885 (2023).
  * Dolgov and Khoromskij (2013) S. Dolgov and B. Khoromskij, SIAM J. Matrix Anal. A. 34, 593 (2013).
  * Durran (2010) D. R. Durran,  _Numerical Methods for Fluid Dynamics: With Applications to Geophysics_ (Springer, 2010).
  * Lasser and Su (2022) C. Lasser and C. Su, J. Math. Phys. 63, 072107 (2022).
  * Dolgov and Savostyanov (2014) S. V. Dolgov and D. V. Savostyanov, SIAM J. Sci. Comput 36, A2248 (2014).
  * Einkemmer and Jospeh (2021) L. Einkemmer and I. Jospeh, J. Comput. Phys. 443, 110495 (2021).
  * Oseledets and Tyrtyshnikov (2010) I. V. Oseledets and E. Tyrtyshnikov, Linear Algebra Appl. 432, 70 (2010).
  * tn_(2023) “Resources for tensor network algorithms, theory, and software,”  (2023).
  * McCulloch (2007) I. P. McCulloch, J. Stat. Mech. , P10014 (2007).
  * Haegeman  _et al._ (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
  * Khoromskij  _et al._ (2012) B. N. Khoromskij, I. V. Oseledets, and R. Schneider, “Efficient time-stepping scheme for dynamics on tt-manifolds,”  (2012), MIS MPI Leipzig:24/2012 .
  * Cappelli  _et al._ (2023) L. Cappelli, F. Tacchino, G. Murante, S. Borgani, and I. Tavernelli, “From Vlasov-Poisson to Schrödinger-Poisson: dark matter simulation with a quantum variational time evolution algorithm,”  (2023), arXiv:2307.06032 [quant-ph] .
  * Frenkel (1934) J. Frenkel,  _Wave Mechanics: Advanced General Theory_ (Clarendon Press, Oxford, 1934) pp. 247–258.
  * Paeckel  _et al._ (2019) S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Ann. Phys. 411, 167998 (2019).
  * Courant  _et al._ (1967) R. Courant, K. Friedrichs, and H. Lewy, IBM Journal of Research and Development 11, 215 (1967).
  * Oseledets and Dolgov (2012) I. V. Oseledets and S. V. Dolgov, SIAM J. Sci. Comput. 34, A2718 (2012).
  * Sonneveld (1989) P. Sonneveld, SIAM J. Sci. Stat. Comp. 10, 36 (1989).
  * Gray (2018) J. Gray, Journal of Open Source Software 3, 819 (2018).
  * Orszag and Tang (1979) S. A. Orszag and C.-M. Tang, J. Fluid Mech. 90, 129–143 (1979).
  * Juno  _et al._ (2018) J. Juno, A. Hakim, J. TenBarge, E. Shi, and W. Dorland, J. Comput. Phys. 353, 110 (2018).
  * Birn and Hesse (2001) J. Birn and M. Hesse, J. Geophys. Res. Space 106, 3737 (2001).
  * Schmitz and Grauer (2006) H. Schmitz and R. Grauer, Phys. Plasmas 13, 092309 (2006).
  * Harris (1962) E. G. Harris, Il Nuovo Cimento (1955-1965) 23, 115 (1962).
  * Yang and White (2020) M. Yang and S. R. White, Phys. Rev. B 102, 094315 (2020).
  * Einkemmer and Lubich (2019) L. Einkemmer and C. Lubich, SIAM J. Sci. Comput. 41, B1061 (2019).
  * Einkemmer  _et al._ (2020) L. Einkemmer, A. Ostermann, and C. Piazzola, J. Comput. Phys. 403, 109063 (2020).
  * Hakim  _et al._ (2020) A. Hakim, M. Francisquez, J. Juno, and G. W. Hammett, J. Plasma Phys. 86, 905860403 (2020).
  * Hakim and Juno (2020) A. Hakim and J. Juno, in  _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ (IEEE Press, 2020).
  * Kazeev  _et al._ (2013b) V. A. Kazeev, B. N. Khoromskij, and E. E. Tyrtyshnikov, SIAM J. Sci. Comput. 35, A1511 (2013b).
  * Kazeev  _et al._ (2017) V. Kazeev, I. Oseledets, M. Rakhuba, and C. Schwab, Adv. Comput. Math. 43, 411 (2017).
  * Kazeev and Schwab (2018) V. Kazeev and C. Schwab, Numerische Mathematik 138, 133 (2018).
  * Kazeev  _et al._ (2022) V. Kazeev, I. Oseledets, M. Rakhuba, and C. Schwab, Multiscale Model. Sim. 20, 893 (2022).
  * Oseledets and Tyrtyshnikov (2011) I. V. Oseledets and E. E. Tyrtyshnikov, SIAM J. Sci. Comput. 33, 1315 (2011), publisher: Society for Industrial and Applied Mathematics.

See pages 1 of QTT_for_Vlasov_SI.pdf

See pages 2 of QTT_for_Vlasov_SI.pdf

See pages 3 of QTT_for_Vlasov_SI.pdf

See pages 4 of QTT_for_Vlasov_SI.pdf

See pages 5 of QTT_for_Vlasov_SI.pdf

See pages 6 of QTT_for_Vlasov_SI.pdf

See pages 7 of QTT_for_Vlasov_SI.pdf

See pages 8 of QTT_for_Vlasov_SI.pdf

See pages 9 of QTT_for_Vlasov_SI.pdf

See pages 10 of QTT_for_Vlasov_SI.pdf

See pages 11 of QTT_for_Vlasov_SI.pdf

See pages 12 of QTT_for_Vlasov_SI.pdf

See pages 13 of QTT_for_Vlasov_SI.pdf

See pages 14 of QTT_for_Vlasov_SI.pdf

See pages 15 of QTT_for_Vlasov_SI.pdf

See pages 16 of QTT_for_Vlasov_SI.pdf

See pages 17 of QTT_for_Vlasov_SI.pdf

See pages 18 of QTT_for_Vlasov_SI.pdf

See pages 19 of QTT_for_Vlasov_SI.pdf

See pages 20 of QTT_for_Vlasov_SI.pdf

See pages 21 of QTT_for_Vlasov_SI.pdf