Dynamical Tensor Train Approximation for Wigner--Fokker--Planck Equation with Applications to Open Quantum Systems

1 Introduction

• WFP computation combines a nonlocal Wigner pseudo-differential operator, naturally treated in Fourier space, with local transport, friction, and diffusion terms.This structural difference makes operator splitting a standard tool for Wigner-type problems [70].Source: §1 · pp. 1–2
• In three spatial dimensions, 64 grid points per position and momentum direction yield about 6.9 × 10^10 unknowns, requiring roughly 550 GB in double precision.Full Wigner-function resolution is useful for phase-space structure, reduced observables, and system–environment interplay.Source: §1 · p. 2
• Dynamical low-rank methods evolve solutions on a prescribed low-rank manifold, whereas step-and-truncate methods compress after a time step and can incur costly intermediate rank growth.Projector splitting breaks the tangent-space projection into simpler substeps.Source: §1 · p. 2
• Prior low-rank work addressed Wigner-type coherent dynamics or classical Fokker–Planck dynamics, but neither line directly treated their interaction in the open quantum WFP setting.The paper identifies this gap as motivation for investigating low-rank structure in WFP.Source: §1 · pp. 2–3
• The method evolves the Wigner function directly in tensor-train form with projector splitting, and treats the nonlocal operator using low-rank potential differences and matrix-free exponential actions.The method is intended for high-dimensional open quantum dynamics and avoids full tensor-product computations.Source: §1 · p. 3

2 Wigner–Fokker–Planck equation

• The WFP equation combines free transport and a potential-associated Wigner pseudo-differential operator with momentum friction, momentum diffusion, cross-diffusion, and spatial diffusion.The Wigner function is real-valued and depends on time, position, and momentum.Source: §2 · p. 3
• With the friction and diffusion terms included, WFP describes an open quantum system; without them, the equation reduces to the Wigner equation for a closed system.For nonnegative Dpp and Dxx, complete positivity is ensured when the coefficients satisfy the Lindblad condition.Source: §2 · pp. 3–4
• In the momentum Fourier representation, the Wigner pseudo-differential term is evaluated by multiplying the transformed Wigner function by the potential difference, then applying an inverse momentum transform.The formal Moyal expansion’s leading term is the classical force term, recovered as ℏ→0.Source: §2 · p. 4

3.1 Tensor train representation of the Wigner function

• Orthogonal TT cores provide canonical forms used by the projector-splitting method, and their bases can be obtained by SVD or QR decompositions.The method uses a prescribed TT rank and allows complex-valued TT entries because the Wigner operator involves complex arithmetic.Source: §3.1 · p. 5
• The tensor-train representation stores the Wigner function with memory cost O(d r^2(Nx + Np)), rather than the full tensor-product cost, when its intrinsic ranks remain small.The variables are interleaved by position and momentum, and the method evolves the solution directly in TT form.Source: §3.1 · pp. 4, 6

3.2 Low-rank representations of the potential difference

• The potential difference is assumed to admit a low-rank tensor representation; the paper develops CP and TT forms, with Tucker also described as a possible alternative.The potential’s TT ranks are independent of the prescribed ranks used for the Wigner function.Source: §3.2 · pp. 6–7
• The equation is split into local operators evaluated in phase space and a global Wigner operator evaluated naturally in Fourier space.Directly applying the nonlocal operator can substantially increase TT ranks, so projector splitting projects its action onto the prescribed low-rank manifold before forming that high-rank intermediate.Source: §3.2 · p. 7

4 Projector-splitting method

• A time step uses 2d center-core steps and 2d−1 bond steps in a left-to-right sweep, totaling 4d−1 local evolution problems.For d = 1, the method coincides with the K-, S-, and L-substep projector-splitting integrator.Source: §4 · pp. 7–8, §4.4 · p. 14

4.1 Spatial center-core step

• The spatial center-core projection yields the same local evolution form for CP and TT potential representations, with distinct environment-index structures.CP uses a shared separable index; TT connects left and right environments through neighboring potential TT indices.Source: §4.1 · pp. 9–10

4.2 Momentum center-core step

• The momentum center-core step has the same local matrix structure as the spatial step, but evolves the L2 contribution in Fourier space.The remaining Fourier-variable dependence is carried by the potential factor and center core.Source: §4.2 · p. 11

4.3 Bond step

• Bond steps evolve a bond matrix using projected L1 and L2 contributions, with the TT potential bond index connecting the environments across the cut.The bond matrix is independent of physical variables, so its projection uses inner products over all phase-space variables.Source: §4.3 · pp. 12–13

4.4 First-order time discretization

• The first-order time discretization applies a matrix-exponential update to projected L2 and an explicit forward Euler increment to projected L1.Forward Euler is consistent with the first-order projector-splitting scheme; applying it directly to L2 is unstable because its eigenvalues lie on the imaginary axis.Source: §4.4 · p. 13

4.5 Matrix-free implementation

• The matrix-free exponential action reduces rank dependence from dense O(r^6) exponential cost to O(nmvϱr^3) when nmv remains moderate, while avoiding O(r^4) matrix storage.For TT potential representations, operator-application costs depend on neighboring potential ranks instead.Source: §4.5 · p. 14

• The spatial and momentum domains are truncated under the assumption that the Wigner function is negligible outside the chosen phase-space region over the time interval of interest.The method uses the same interval in each spatial direction and another shared interval in each momentum direction for simplicity, though this is not required.Source: §5 · p. 15
• The discrete Wigner function is a 2d-way tensor with alternating spatial and momentum modes, approximated in TT form to reduce memory use.The prescribed solution TT ranks are fixed in this discretization.Source: §5 · p. 15
• Conservative flux discretization keeps transport, diffusion, and drift local to their physical variables, allowing their action to be incorporated directly into the corresponding TT cores.The transport term uses second-order upwind reconstruction in the interior, with zero-flux boundary conditions.Source: §5 · p. 15
• The momentum grid determines the conjugate Fourier grid, and FFTs with a phase correction account for the cell-centered momentum nodes.The same Fourier grid is used for the momentum dependence of the potential difference.Source: §5 · pp. 16–17

6.1 Computational cost for d ≥2

• For fixed d = 2 or 3, one time step scales linearly with Nx and Np, up to the FFT factor log Np, and cubically with solution rank r when nmv and ϱ remain moderate.The local differential-operator contractions scale as O(dr^3(Nx + Np)) and do not involve the CP rank ϱ.Source: §6.1 · p. 18
• For d ≥ 2, CP and TT potential representations have the same dependence on solution rank and grid resolution, while their potential-representation costs depend respectively on ϱ and neighboring TT ranks bounded by λ^2.The matrix-free CP environment construction has cost O(dϱr^3(Nx+Np)+dϱr^2Np log Np).Source: §6.1 · pp. 17, 19

6.2 Computational cost for d = 1

• For d = 1, the TT potential representation has the same complexity estimate as the CP representation after replacing ϱ with the single potential bond rank λ1.The one-dimensional center-core exponential operators have dimension r × r rather than r^2 × r^2, making this case substantially cheaper than d = 2 or 3.Source: §6.2 · p. 19

7 Numerical examples

• The time step in the experiments followed a CFL-type stability bound for explicitly treated L1; the advanced exponential treatment of L2 added no explicit stability restriction.The bound is described as a practical forward-Euler stability estimate, not a sharp condition for the full projector-splitting scheme.Source: §7 · p. 19

7.1 1D1P Wigner–Fokker–Planck equation

• The one-dimensional phase-space results were consistent with [70], validating the method with coherent Wigner transport and environmental friction and diffusion present simultaneously.The comparison showed the solution at t = 0, 2, 4, and 6.Source: §7.1 · p. 20

7.2 2D2P Wigner equation

• Across prescribed TT ranks, phase-space plots were visually similar and differences from the highest-rank reference decreased as rank increased; rank (30, 30, 30) produced evolution qualitatively consistent with [54].The rank-(30, 30, 30) solution was used for the longer-time simulation up to t = 6.Source: Figure 2 · p. 23, §7.2 · p. 23, Table 2 · p. 24

7.3 2D2P open-system dynamics with a weakly nonquadratic potential

• With increasing friction and diffusion, the two initially separated wave packets were progressively damped and merged around the low-potential region; the potential difference had an exact six-term CP representation.The stronger damping occurred in case (c), which used larger friction and diffusion coefficients than case (b).Source: §7.3 · p. 25, Figure 4 · p. 26
• Introducing Fokker–Planck terms substantially reduced the numerical effective TT ranks, consistent with increased compressibility in this example but not establishing a general lower-rank property for open dynamics.The reported interpretation is that dissipation and diffusion suppress fine, oscillatory phase-space structures; the authors limit the conclusion to some regimes.Source: §7.3 · pp. 25, 27, Figure 5 · p. 27

7.4 3D3P simulation of an open quantum system

• In the six-dimensional open-system test, the tensor-train representation reduced storage requirements by approximately five orders of magnitude, demonstrating feasibility in this high-dimensional regime.The discretization used 64 grid points in each of six directions and prescribed TT rank (24, 24, 24, 24, 24).Source: §7.4 · p. 28
• The reduced Wigner function and spatial density evolve similarly under the two CFL coefficients, indicating consistency with time-step refinement.These reduced quantities can be evaluated from the TT representation without reconstructing the full six-dimensional tensor.Source: §7.4 · p. 28, Figure 7 · p. 29
• The real mass remains close to unity and the imaginary mass close to zero; the smaller CFL coefficient produces a smaller mass error.Although the equation conserves total mass analytically, the numerical method does not impose exact mass conservation.Source: §7.4 · p. 29, Figure 8 · p. 29
• Unlike density-matrix low-rank schemes that enforce complete positivity and trace preservation [22], this method imposes neither property discretely; conservative tensor-train formulations remain future work.The small mass errors in this experiment nevertheless indicate accurate maintenance of trace normalization.Source: §7.4 · p. 29
• Runtime scales empirically as O(N^1.07r^2.12), consistent with the complexity estimate and FFT-related costs over the tested grid and rank range.Timings average one propagation step over 20 consecutive steps; the nearly linear grid and approximately quadratic rank dependence are associated with potential-environment FFT operations.Source: §7.4 · pp. 29–30, Table 3 · p. 30

8 Conclusion

• The method evolves the Wigner function directly in tensor-train form with projector splitting, and treats the nonlocal operator using low-rank potential-difference representations and matrix-free exponential actions.Other Wigner–Fokker–Planck terms are evaluated and projected directly in phase space, avoiding construction of the full phase-space tensor.Source: §8 · p. 30
• Numerical experiments in one-, two-, and three-dimensional settings demonstrate the method’s accuracy and efficiency.The three-dimensional setting corresponds to a six-dimensional position–momentum phase space.Source: §8 · p. 30
• The current method uses fixed tensor-train ranks and a first-order projector-splitting time discretization, and does not enforce trace preservation or complete positivity at the discrete level.Future work includes rank-adaptive and higher-order methods, structure-preserving tensor trains, and comparisons with alternative tensor formats and integrators.Source: §8 · p. 30