Tensor Network Compression for Fully Spectral Vlasov-Poisson Simulation

TL;DR

The paper tackles the challenge of efficiently simulating collisionless kinetic plasmas by representing the six-dimensional phase-space distribution $f(x,v,t)$ in a quantics tensor-train (TT) format and evolving it with a fully spectral, compressed solver. Time stepping uses a Strang splitting that alternates advection and acceleration substeps, with spectral transforms applied directly in TT form via TT–MPOs and tensor cross interpolation (TCI) for the nonlinear substep, while the self-consistent electric field is computed within the TT framework. The authors demonstrate the method on 1D1V Landau damping and the two-stream instability, showing accurate damping/growth rates and robust conservation properties, while systematically exploring how TT truncation and bond dimensions affect fidelity, negativity artifacts, and runtime. The results indicate that tensor-network compression can provide substantial memory and computational advantages for high-fidelity Eulerian Vlasov simulations, with clear directions for improving positivity, extending to higher dimensions, and optimizing performance.

Abstract

We propose a numerical method for kinetic plasma simulation in which the phase-space distribution function is represented by a low-rank tensor network with an adaptive level of compression. The Vlasov-Poisson system is advanced using Strang splitting, and each substep is treated spectrally in the corresponding variable. By expressing both the distribution function and the Fourier transform as tensor network objects (state and operator representations), spectral transforms are applied directly in compressed form, enabling time stepping without reconstructing the full phase-space grid. The self-consistent electric field is also computed within the tensor formalism. The charge density is obtained by contracting over velocity degrees of freedom and extracting the zero Fourier mode, which provides the source term for a spectral Poisson solver. We validate the approach on standard benchmarks, including Landau damping and the two-stream instability. Finally, we systematically study how compression parameters, including truncation tolerances and internal ranks (bond dimensions), affect momentum and energy conservation, positivity behavior, robustness to filamentation, and computational cost.

Figures (9)

• Figure 1: Workflow for the spectral Strang splitting integration scheme described in Eq. \ref{['eq:strang_splitting_integration']} over $N_s$ time steps. For both advection and acceleration, the solution is transformed to the Fourier domain of the variable of differentiation, after which the corresponding time propagator is applied. The highlighted block indicates the repeated inner loop applied $N_s-1$ times.
• Figure 2: (a) The discretized distribution function $f(x_{i_x}, v_{i_v})$ is encoded as a rank-$2R$ tensor using the quantics representation and subsequently compressed into a tensor train (TT) format. In a tensor network diagram, nodes denote tensors, while edges connecting two nodes indicate summation (contraction) over a shared index. Edges that terminate without connecting to another node represent free (site) indices of the tensor. (b) Application of a matrix product operator (MPO) to the TT representation of $f$, yielding a new TT with the same physical structure. The action of the MPO is performed by contracting the site indices of the TT with the corresponding input indices of the MPO at each site.
• Figure 3: Time evolution of the TT representation of the initial condition $f^0$, following the sequence dictated by the second-order Strang splitting scheme. Advection steps are implemented via TT--MPO contractions, whereas the acceleration steps are realized using tensor cross interpolation (TCI). Fourier transforms in position (red) and velocity (yellow) act only on the corresponding TT site registers, whereas the phase shifts associated with advection and acceleration (green and blue) act on the full position--velocity register.
• Figure 4: (a) The electron density $\rho_e$ is computed by contracting the tensors in the velocity register by $10^T$, effectively selecting the zero-mode. One site must be contracted by $\sqrt{2^R}\,\Delta v0^T$, to consider the spacing of the spatial grid and to compensate for the prefactor of the unitary Fourier transform. (b) The electric field is acquired by multiplying the TT representation of the charge density $\rho$ by a sequence of three MPOs, corresponding to Eq. \ref{['eq:poisson_fourier_kernel']}.
• Figure 5: TT-based simulation of Landau damping, with initial condition set according to Eq. \ref{['eq:ic_landau_damping']}. (a) Energy of first mode of electric field $\mathcal{E}_\text{field}^{(1)}$. The analytically computed damping rate is indicated by the dashed line, showing good agreement with the simulated results. (b--c) Evolution of absolute momentum $\lvert P \rvert$ and absolute relative energy deviation $\lvert \mathcal{E}(t) -\mathcal{E}(0) \rvert / \lvert \mathcal{E}(0) \rvert$ throughout the simulation. Whereas the momentum appears to fluctuate, the energy conservation quickly stabilizes. (d--f) Snapshots of the distribution function $f$, showing gradual filamentation along the velocity domain.