Pseudospectral method for solving PDEs using Matrix Product States
Jorge Gidi Paula García-Molina Luca Tagliacozzo Juan José García-Ripoll







jorgegidi@udec.cl
Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b, Madrid 28006, Spain
Millennium Institute for Research in Optics and Departamento de Física, Facultad de Ciencias Físicas y Matemáticas, Universidad de Concepción, Casilla 160-C, Concepción, Chile


paula.garcia@iff.csic.es
Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b, Madrid 28006, Spain

Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b, Madrid 28006, Spain

Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b, Madrid 28006, Spain


§ ABSTRACT 
  
  This research focuses on solving time-dependent partial differential equations (PDEs), in particular the time-dependent Schrödinger equation, using matrix product states (MPS). We propose an extension of Hermite Distributed Approximating Functionals (HDAF) to MPS, a highly accurate pseudospectral method for approximating functions of derivatives. Integrating HDAF into an MPS finite precision algebra, we test four types of quantum-inspired algorithms for time evolution: explicit Runge-Kutta methods, Crank-Nicolson method, explicitly restarted Arnoli iteration and split-step.
  The benchmark problem is the expansion of a particle in a quantum quench, characterized by a rapid increase in space requirements, where HDAF surpasses traditional finite difference methods in accuracy with a comparable cost.
  Moreover, the efficient HDAF approximation to the free propagator avoids the need for Fourier transforms in split-step methods, significantly enhancing their performance with an improved balance in cost and accuracy. 
  Both approaches exhibit similar error scaling and run times compared to FFT vector methods; however, MPS offer an exponential advantage in memory, overcoming vector limitations to enable larger discretizations and expansions. Finally, the MPS HDAF split-step method successfully reproduces the physical behavior of a particle expansion in a double-well potential, demonstrating viability for actual research scenarios.



Pseudospectral method for solving PDEs using Matrix Product States
    Juan José García-Ripoll
    Received ????; accepted ????
==================================================================




  

§ PRIVATE: MAIN MESSAGES – TODO RM

  What messages do we plan to deliver?
  
    
  * MPS can be more efficient than standard arrays to find the time evolution in practical, large physical problems
    
  * We propose the HDAF as an efficient and precise tool to approximate functions of derivatives in the context of MPS-MPO
    
  * The free evolution operator is efficiently represented in the coordinate basis using the HDAF-MPO formalism
    
  * We can implement split-step evolution without Fourier transforms and achieve subexponential time scaling.
    
  * We demonstrate the convergence of the method competitive to vectors in particle expansion problems.
    
  * In the absence of chirping, our results are favorable to MPS in contrast to vectors using FFT.
  




§ INTRODUCTION


Solving time-dependent partial differential equations (PDEs) is crucial across most fields in science and engineering.
In the quantum domain, the challenge is additionally plagued by exponential costs <cit.>, both in the number of components of the system and in the unbounded domain of certain problems. An example of the latter is a particle expansion in a potential well <cit.>.


The simulation of large quantum systems rapidly becomes untractable by traditional computational approaches due to their exponential scaling in complexity. Quantum computers have been posed as a promising tool to solve PDEs more efficiently <cit.>, making use of the exponential compression that the amplitude encoding provides <cit.> and tools such as the Quantum Fourier Transform (QFT) <cit.>. However, there are still no scalable and fault-tolerant quantum computers where those algorithms may be applied <cit.>, and it is also uncertain if quantum computers are required, especially for bandwidth-limited functions with limited entanglement <cit.>.


A current challenge is to develop alternative, quantum-inspired algorithms that profit from some of the exponential compressions available in quantum algorithms and encodings while enabling the execution in classical computers.
This challenge involves developing three tools: (i) an efficient encoding of the PDE solution, (ii) a corresponding encoding of the PDE operator itself, and (iii) an algorithm for time evolution.
In this work, we address these three objectives in novel ways.
We use the representation of matrix product states (MPS) for bandwidth-limited functions <cit.>, also known in the mathematical community as quantized tensor trains (QTT) <cit.>. Within this formulation, we develop an innovative encoding of differential operators, PDEs, and free-evolution propagators based on Hermite Distributed Approximating Functionals (HDAF) <cit.>, to yield accurate matrix product operators (MPOs) with a low bond dimension. We present three families of time evolution algorithms using these operators: Global explicit and implicit methods, Arnoldi-based methods, and split-step methods. These new techniques are benchmarked against each other and alternative quantum-inspired methods based on finite differences <cit.>. They are also compared to state-of-the-art spectral split-step methods in the standard vector representation <cit.>.


As a benchmark, we study the expansion of a particle in a potential well. This computationally demanding scenario is highly relevant in optomechanics research <cit.>. As the domain size increases dramatically, conventional computational methods suffer considerable strain. In this context, the qubit encoding in MPS/QTT may provide exponential data compression. However, the particle's acceleration induces chirping of the wavefunction, potentially increasing the bond dimension required beyond any computational benefit. For the limited cases where this problem can be analytically solved <cit.>, it constitutes an effective testbed to stress and validate numerical solvers.


When correctly tuned, the HDAF MPOs for differential operators are much more accurate than those based on finite differences and have a comparable cost. Also, the time evolution methods built upon these operators inherit an improvement in accuracy and efficiency.
Our best-performing time evolution method is the HDAF split-step, which leverages the efficient HDAF representation of the free propagator on a coordinate basis to avoid using Fourier Transforms.

For the particle expansion problem, quantum-inspired methods converge with subexponential time scaling, which is competitive with a vector implementation. Moreover, in the absence of chirping, our results are favorable to the HDAF split-step method using MPS, in contrast to vectors using the Fast Fourier Transform (FFT).

Quantum-inspired methods have been used to tackle quantum numerical analysis problems before, such as static <cit.> and time-dependent PDEs <cit.>, and they have even permeated to other fields, such as kinetic plasma simulation <cit.>. Time evolution problems usually rely on Fourier techniques with Trotter expansions <cit.> or Chebyshev propagation schemes <cit.>. As an alternative, this work extends the HDAF approach to MPS/QTT, accurately and efficiently encoding differential operators, such as arbitrary functions of derivatives. Applied to the free-propagator, it enables Trotter expansions without Fourier transforms.


The work is structured as follows.
Section <ref> presents the expansion problems to validate the numerical methods.
Section <ref> introduces the MPS finite-precision algebra framework for quantum-inspired numerical algorithms and the HDAF machinery to approximate derivatives and functions of derivatives as MPOs, including metaheuristics to tune these approximations.
Section <ref> reviews the quantum time evolution problem and presents the numerical methods to address it. These techniques are contrasted via a one-step study on the one-dimensional quantum quench problem to determine the most convenient option.
Then, Section <ref> examines the best-performing technique for a long-time evolution on the complete range of expansion for the harmonic and non-harmonic problems.
Finally, Section <ref> summarizes the conclusions of this study.




§ PARTICLE EXPANSION


A quantum quench is a fundamental process where a system is driven out of equilibrium by a sudden change in its Hamiltonian <cit.>. In this section, we introduce the problem of particle expansion due to the sudden relaxation of an initial harmonic potential.

The particle expansion process presents several characteristics relevant to our study. First, it is a problem of interest in many areas, including many-body physics <cit.> and quantum optomechanics <cit.>. Second, it stresses traditional and MPS-based numerical methods. Vector representations become too expensive for large expansions as they require the storage of a huge spatial domain. This makes the problem an interesting playground for MPS simulation. Nonetheless, the particle acceleration induces a chirping of the wavefunction, thus incrementing the bond dimension of the solution and putting a strain on MPS-based simulations as well. A third feature of this process is that when the quench varies from one harmonic potential to another, a closed-form analytical solution is known <cit.>. This is a critical feature allowing us to benchmark not only the speed but also the accuracy of our methods.


The Schrödinger equation describes the evolution of a particle's wavefunction

    i∂_t ψ(x,t) = (-ħ^2/2m∂_x^2 + V(x,t))ψ(x,t).


The election of the potential V(x,t) determines the physical behavior of the system. The benchmark problem we introduce is the sudden change of a harmonic potential from frequency ω_0 to ω_H at time t=0,




    V(x,t) = 
      1/2ω_0^2x^2,    t ≤ 0,
    1/2ω_H^2x^2,    t> 0.


Assuming the wavefunction at time t = 0 is relaxed to the ground state of the previous Hamiltonian,

    ψ(x, t=0)    = (ω_0/π)^1/4exp(1/2ω_0x^2),

its evolution is prescribed according to

    ψ(x,t)    = (ω(t)/π)^1/4exp(-[ω(t)/2 + iβ(t)] x^2), 
    ω(t)    = ω_H(ω_H/ω_0cos^2(ω_H t)+ω_0/ω_Hsin^2(ω_H t))^-1,
    β(t)    = ω(t)/4(ω_H/ω_0-ω_0/ω_H)sin(2ω_H t).

The solution is a complex Gaussian with width σ(t)=1/√(ω(t)), modulated in time with period π/ω_H. For ω_H < ω_0, the system undergoes an expansion during the first half of the period. The smallest width is σ_min = 1 / √(ω_0) at time t=0, and the largest width σ_max = √(ω_0) / ω_H at time t = 0.5π / ω_H.

The expansion is quantified by σ_max / σ_min = ω_0 / ω_H. That is, the frequency ratio ω_0/ω_H is also the expansion ratio, dictating the total amplification of the wavefunction's spatial extent.
A larger expansion ratio requires more points to represent the solution accurately, challenging common computational approaches.

Another problem featuring similar expansion dynamics occurs when the initially harmonic potential changes to a wider trap with a non-harmonic component. This scenario holds high interest for experimental settings in optomechanics research <cit.>. However,
analytical solutions are not usually known, while the intricacies that complicate numerical methods still hold.
One such case happens for a double-well potential,

    V(x,t) =
      1/2ω_0^2x^2,     t ≤ 0,
    1/2ω_H^2 x^2 + u exp(-x^2/2σ^2),     t> 0.


Qualitatively, the harmonic expansion behavior dominates the wavefunction's evolution for a sufficiently small value of u, with the Gaussian bump at x=0 acting as a perturbation. This allows us to use the results of the harmonic quench analysis as a basis for the study of the double-well potential. The Gaussian component of the potential, with u > 0, is expected to induce a symmetric separation of the expanding particle. This will deviate the wavefunction's evolution from the Gaussian form (<ref>) into a state with a two-peaked probability density. We choose this problem as a test to account for the feasibility of using the proposed numerical methods against a setting of actual research interest.




§ MPS ENCODING AND DIFFERENTIAL OPERATORS


MPS originally arose in the domain of physics to study quantum many-body problems. Dolgov rediscovered this formalism in the field of applied mathematics under the name of quantized tensor trains (QTTs) <cit.>, a subclass of the broader class of tensor trains (TTs) <cit.>. A continuous alternative class of TTs is the functional tensor trains (FTTs) <cit.>, which hold a network of univariate matrix-valued functions instead of rank-three tensors. 

The exponential memory compression of MPS/QTT can bypass the curse of dimensionality in the representation of functions. This motivated the development of similar encodings from a quantum-inspired perspective <cit.>. Indeed, MPS/QTT constitute efficient ansätze for representing functions with fastly decaying Fourier coefficients <cit.>.

MPS and TTs have been successfully applied to a variety of numerical analysis problems such as high-dimensional nonlinear PDEs <cit.>, the Hamilton Jacobi Bellman equations <cit.>, the Schrödinger equation <cit.>, and stochastic problems <cit.>. Combining Fourier techniques with Trotter expansions <cit.>, or Chebyshev propagation schemes <cit.> allows for quantum dynamics simulations. Other approaches rely on one-step implicit time integration using an ALS-type solver or global space-time formulation to solve multi-dimensional parabolic problems <cit.>. Quantum-inspired proposals are successful at solving different PDEs, such as Schrödinger equations <cit.>, turbulence problems <cit.>, Hamiltonian PDEs <cit.>, and the Vlasov-Poisson system <cit.>.

Combining this MPS-based representation of functions and the analogous encoding of operators as matrix product operators (MPO), along with a basic set of operations and efficient truncation algorithms, leads to a finite precision algebra <cit.> that operates similarly to standard matrix-vector operations. This algebra is the basis for developing time evolution quantum-inspired algorithms in section <ref>.

This section presents an extension of the Hermite Distributed Approximating Functionals (HDAF) to reconstruct functions of differential operators within this finite precision MPS-MPO framework. The HDAF performs this reconstruction as a linear combination of Hermite polynomials weighted by a Gaussian filter, approximating these operators with tunable pseudospectral precision at a limited cost. The section also covers a general review of the HDAF and the metaheuristics behind the use of this technique.



 §.§ Quantum-inspired numerical analysis


The amplitude encoding of functions <cit.> represents a function f(x) with x∈ [a,b) in an n-qubit quantum register as a normalized quantum state

    |f^(n)⟩ = 1/𝒩_f^1/2∑_s=0^2^n-1 f(x_s^(n))|s⟩,    x_s^(n) = a + sΔ x ^(n),

where 𝒩_f is the normalization constant and the index s={ 0,…, 2^n} labels each coordinate x_s to its corresponding quantum state |s⟩ <cit.>. The amplitude encoding may be exponentially compressed in a quantum register by mapping the indices s to the states of n qubits. This leads to a binary encoding of the coordinates s=∑_i=1^n2^n-is_i, with s_i = { 0, 1 }.

This function representation in a quantum register creates a many-body wavefunction that, for bandwidth-limited functions, admits an efficient MPS representation <cit.>

    |ψ⟩   =∑_{ s_i }∑_{α_i } (A_α_1^s_1A_α_1,α_2^s_2… A_α_N-1^s_N)|s_1⟩
                 ⊗|s_2⟩⊗ ... ⊗|s_N⟩,

where each qubit is mapped to a site of the MPS. The physical indices s_i correspond to the values of the energy levels of the qubits, and α_i are the bond indices, where α_i = {1,…,χ_i} with χ_i the bond dimension.


Given this efficient representation of functions, a relevant challenge is to find a comparable representation for operators that act on these functions, from PDE operators that describe the action of derivatives and potentials to evolution operators that model the dynamics of the quantum state.
Potentials are diagonal operators and can be derived from the MPS representation of their corresponding functions either exactly or via interpolation techniques like the Chebyshev approximation <cit.> or the TT-cross interpolation <cit.>. Thus, the real challenge arises from
finding suitable, efficient representations of derivatives and functions thereof. 

The finite difference method is one of the most widespread techniques for approximating derivatives. This method uses the order p Taylor expansion of the function, with an error O(Δ x^m), where m depends on the terms combined for the approximation. Thus, the grid size limits the accuracy of the finite difference method. The most common implementations are the centered finite difference formulas

    ∂ f(x)/∂ x   = f(x + Δ x) - f(x - Δ x)/2 Δ x + O(Δ x^2), 
    ∂^2 f(x)/∂ x^2   = f(x + Δ x) - 2 f(x) + f(x - Δ x)/Δ x^2 + O(Δ x^2),

whose truncation error scales quadratically with the discretization. The weak noise suppression of the centered finite difference formula can be enhanced to construct smooth noise-robust differentiators <cit.>.

The finite difference MPO uses a linear combination of displacement operators Σ^±<cit.>,

    |∂_xf^(n)⟩   ≃1/2Δx(Σ̂^+-Σ̂^-)|f^(n)⟩ + O(Δ x^2), 
    |∂^2_xf^(n)⟩   ≃1/Δx^2(Σ̂^+-2𝕀+Σ̂^-)|f^(n)⟩ + O(Δ x^2).

This MPO representation of the finite difference operators has a fixed bond dimension χ=3 independent of the number of sites.

In addition to the truncation error, the round-off error affects the finite-difference scheme. The trade-off of truncation and round-off error determines the optimum step size Δ x since the truncation error decreases with Δ x, but still, values that are too small will result in numerical errors due to round-off. Let us consider the second-order derivative approximation since this appears in the propagator of the time evolution. 
The round-off error occurs when the difference between two consecutive points |f(x)-f(x±Δ x)| of the discretized function is of the order of the machine precision δ. This error is proportional to δ/Δ x^2 and can be amplified by the small denominator Δ x^2. Increasing the step size while maintaining the number of points for the discretization can correct this error. 





 §.§ Hermite Distributed Approximating Functionals


This work aims to overcome finite difference formulas' precision and speed limitations. In this regard, it is well known that spectral methods, such as Fourier techniques, can provide exponential speedup and precision guarantees for sufficiently smooth functions. Such methods have been derived in the domain of MPS/QTT methods <cit.> with some ad-hoc heuristics to improve the creation of operators. Hermite Distributed Approximating Functionals (HDAF) give a powerful yet lesser-known numerical analysis technique. In the following, we will show how these methods work and how they can be adapted to engineer MPOs of arbitrary differential operators with relatively low costs. Our work differs from previous studies where they have been applied in matrix form to each site of a tensor train <cit.>.


The Distributed Approximating Functionals (DAF) are well-tempered approximations to the Dirac delta distribution.
The first DAF, developed before the class name was coined, was the Hermite DAF <cit.>,

    δ_M(x; σ) = exp(-x^2/2σ^2)/√(2π)σ∑_m=0^M/2(-1/4)^mH_2m(x/√(2)σ)/m!,

where H_n(x) is the n-th Hermite polynomial. It has two free parameters: The even integer M and the positive real σ are the order of the highest polynomial and the width of the approximation to the delta distribution, respectively.

The kernel defined in Eq. (<ref>) is a nascent delta function that operates as the identity for polynomial functions of degree M or lower,

    f(x) ≈∫ dx' δ_M(x - x'; σ) f(x'),

approaching the Dirac distribution in the limit σ/M → 0. However, unlike the exact delta distribution, δ_M(x; σ) is generally a bandwidth-limited, infinitely smooth function amenable to quadrature methods and differentiation.

The method's well-tempered property arises from the absence of special points in the reconstruction. Exact reproduction at grid points is not required; therefore, it does not constitute an interpolation scheme.
Moreover, a fundamental property of Eq. (<ref>) is that the approximation converges uniformly to f(x) <cit.>, and the approximation error usually resembles the function f(x) albeit several orders of magnitude smaller <cit.>. This implies that the error is smaller when the function approaches zero, which is relevant and desirable for our wavefunctions applications.

Equation (<ref>) is customarily discretized on a uniform grid with spacing Δ x, using midpoint integration to render it as a matrix-vector product, f(x_i) = ∑_jK_ijf(x_j), where K is a symmetric Toeplitz matrix with components

    K_ij = Δ x δ_M(Δ x |i - j|; σ).


Since δ_M(x; σ) has an exponentially decaying envelope, two essential properties arise: (i) the reconstruction matrix K can be made highly sparse and concentrated around its main diagonal, with the number of diagonals controlled by σ / Δ x, and (ii) a minimal number of diagonals (i.e. quadrature nodes) are required to accurately discretize the integral in Eq. (<ref>). The composite midpoint rule converges especially fast for periodic or peaked functions with vanishing derivatives at the integration limits <cit.>. Moreover, in this context, the midpoint rule surpasses the accuracy of higher-order Newton-Cotes quadrature rules using the same number of grid points <cit.>.


The MPO corresponding to the K matrix can be constructed as a weighted combination of the displacement operators Σ^± as

    K̂   = Δ x δ_M(0; σ) 𝕀
            +  ∑_i=1^2^n - 1Δ x δ_M(iΔ x; σ) (Σ̂^+i + Σ̂^-i),

where the symmetry of δ_M has been used. While in principle the sum ranges over the whole grid, in practice, one only needs to sum until δ_M has vanished according to a prescribed tolerance, as detailed on Sec <ref>.



 §.§ HDAF differentiation



The HDAF formalism opens a path to estimate a function's derivative of any order, as well as functions D[∂/∂ x] of such derivatives. From Eq. (<ref>) it follows that

    D[∂/∂ x] f(x) ≈∫ dx'  D[∂/∂ x] δ_M(x - x') f(x').

An analytical expression for D[∂/∂ x] δ_M(x - x') is usually easy to find. The typical procedure involves using the Rodrigues formula for the Hermite polynomials,

    H_n(x) = (-1)^nexp(x^2)∂^n/∂ x^nexp(-x^2),

to rewrite Eq. (<ref>) as

    δ_M(x; σ)    = 1/√(2π)σ∑_m=0^M/21/m!(-σ^2/4)^m
                   ×∂^2m/∂ x^2mexp(-x^2/2σ^2).

Then, the differential operator is applied over δ_M(x - x'; σ) acting on the exponential, and the Rodrigues formula (<ref>) is used to recover an expression in terms of Hermite polynomials, without explicit derivatives. For instance, the l-th derivative leads to

    δ_M^(l)(x; σ)    = (-1/√(2)σ)^lexp(-x^2/2σ^2)/√(2π)σ
           ×∑_m=0^M/2(-1/4)^mH_2m + l(x/√(2)σ)/m!.


Analog to the reconstruction (<ref>), the differentiating MPO for the l-th derivative in the HDAF formalism is

    K̂^(l)   = Δ x δ_M^(l)(0; σ) 𝕀
          +  ∑_i=1^2^n - 1Δ x δ_M^(l)(iΔ x; σ) (Σ̂^+i + (-1)^lΣ̂^-i),

where the symmetry or antisymmetry of δ_M^(l) has been used for l even or odd, respectively.

Attention must be paid to the fact that the l-th derivative has a prefactor σ^-(l+1) with σ = O(Δ x). Since these operators will be used within a finite-precision numerical framework, round-off errors can significantly impact the accuracy of the differentiation. These errors arise from the limited significant digits that computers may represent, inducing small approximation errors amplified by large weighting prefactors. In particular, these deviations dominate when Δ x is too small, thus limiting how dense the numerical grid can be.

The round-off error amplification problem is expected in numerical differentiation techniques but can be mitigated.
In the particular MPS-MPO framework, differentiating operators can be discretized to accommodate a certain number of qubits and then be adapted to a denser grid, where identities are added as new sites to account for each additional qubit. This approach is equivalent to nearest-neighbor interpolation. It keeps the round-off errors constant and does not introduce extra complexity to the MPO.

Figure <ref> accounts for the differentiation accuracy of HDAF and finite differences in approximating the second derivative of a Gaussian function. In the case of the HDAF, the convergence is faster in the number of qubits for larger values of M, as expected. However, HDAF and finite differences suffer from relevant round-off errors for many qubits. Once the optimal accuracy is achieved for a certain number of qubits, the accuracy of the operators can be retained while acting on a finer grid. Round-off errors are effectively kept constant by following the procedure above.







 §.§ HDAF free propagator


The HDAF scheme is very powerful, as it cannot only approximate derivatives but also general functions of those derivatives. The first and one of the most relevant applications of this idea, posed in Ref. <cit.>, is the
banded approximation of the free propagator.

From Eq. (<ref>), taking the differential operator D[∂/∂ x] to be the free propagator,

    T(τ) = exp(-iτ/2∂^2/∂ x^2),

the kernel for the approximation is T(τ)δ_M(x - x'; σ).
This quantity is readily computed from Eq. (<ref>) since T(τ) commutes with the derivatives and it spreads a Gaussian function to another,

    T(τ)exp(-(x - x')^2/2σ^2)    = [σ/σ_t]exp(-(x-x')^2/2σ_t^2),

mapping the original variance σ^2 to σ_t^2 = σ^2 + iτ. Then, using Eq. (<ref>) yields the free propagator kernel

    δ_M(x - x'; σ, τ)    = T(τ) δ_M(x - x'; σ) 
       = exp(-(x-x')^2/2σ_τ^2)/√(2π)σ_τ
         ×∑_m=0^M/2(-σ^2/4σ_τ^2)^mH_2m(x - x'/√(2)σ_t)/m!.



The MPO for the free propagator in the HDAF formalism is

    K̂_τ   = Δ x δ_M(0; σ, τ) 𝕀
            +  ∑_i=1^2^n - 1Δ x δ_M(iΔ x; σ, τ) (Σ̂^+i + Σ̂^-i).



The kernel (<ref>) becomes complex and highly oscillating as time increases. Also, its width increases as a fundamental consequence of the free propagator. However, the spreading in the HDAF formalism is the minimum possible since it is inherited from the Gaussian generator of the Hermite polynomials <cit.>.

The HDAF approximation for the propagator has been used in many applications of split-step integration methods within the traditional vector framework. In that framework, K̂_τ is represented as a matrix acting on a discretized function. A central contribution in this work is to realize that the same matrix can be more efficiently represented as an MPO using the displacement operators Σ̂^± and additional simplification steps that significantly reduce the effective bond dimension of the operator. In this scenario, the MPO HDAF propagator is a competitive alternative to using MPO Fourier-based techniques <cit.>, directly representing the evolution operator in the coordinate representation.





 §.§ HDAF metaheuristics




  §.§.§ Free parameter election






  
Identical reconstruction.

The formulation of the HDAF operator (<ref>) as a discrete matrix has 2 sources of error: (i) The assumption that the function f(x) can be expressed as a polynomial of degree M under the extent of the Gaussian envelope of the HDAF, and (ii) the discretization of the convolution integral to a finite sum employing the midpoint rule. While the error (i) vanishes in the limit M/σ→∞, it is clear from (ii) that it is not possible to increase M or decrease σ indefinitely. A more oscillatory integrand will require a larger value of σ / Δ x for the Gaussian envelope to cover enough nodes and achieve satisfactory integration accuracy.

In general, for a fixed M, there is a value of σ/Δ x that makes the reconstruction optimal, and the larger is M, the better the maximum accuracy that can be achieved. The rationale behind this optimal relationship between M and σ/Δ x is that a perfect reconstruction happens when the M zeroes of the HDAF match with the zeroes on the grid, and only the term in the origin contributes <cit.>. One possible approach, therefore, is to set the discrete HDAF to be 1 at the origin <cit.>,

    K_ii = Δ x δ_M(0; σ) = 1,

yielding,

    σ_M = Δ x/√(2π)∑_m=0^M/2(-1/4)^mH_2m(0)/m!,

which makes the HDAF approximately vanish at integer multiples of Δ x <cit.>.


In practice, a lower bound to σ/Δ x is prescribed to ensure convergence of the midpoint rule when M is small. For the context of double floating-point precision, we heuristically set

    σ / Δ x ≥ 3
      ⇒σ_min = 3Δ x,

and choose the value of σ according to

    σ = max( σ_M, σ_min).




  
Differentiation.

Assuming that the l-th derivative of the function f(x) is accurately described with an HDAF of order M, i.e., also pertains to the DAF-class <cit.>, then the approximation to the derivative can be thought as a reconstruction of f^(l)(x) instead of f(x),

    f^(l)(x)    ≈∫ dx' δ_M^(l)(x - x'; σ) f(x') 
       = ∫ dx' δ_M(x - x'; σ) f^(l)(x').


In this spirit, the optimal value of σ is again computed using Eq. (<ref>), which does not depend upon the function to reconstruct, provided M is fixed.




  
Free evolution.

Since the action of the free propagator is to spread the original wavefunction, the width of the HDAF will not be a problem for the midpoint integration. For efficiency purposes, the election of σ is made such that the new width of the freely-propagated HDAF is the smallest possible <cit.>. This value follows from equation (<ref>). The leading Gaussian,

    exp(-x^2/2(σ^2 + iτ)) =
      exp(-x^2/2 w^2)
      exp(ix^2/2w^2τ/σ^2),

has an effective variance w^2 = (σ^2 + τ^2/σ^2) with an optimal value σ = √(τ) that minimizes its spatial extent. Then, the value of σ is chosen

    σ = max( σ_M, σ_min, √(τ)),

where σ_M and σ_min and are the same values (<ref>) and (<ref>) used for identical HDAF reconstruction.




  §.§.§ Self-consistent error estimation



The HDAF filter (<ref>) can be analyzed in Fourier space. From equation (<ref>), the kernel spectrum has the analytical form

    δ_M(k; σ) = exp(-k^2σ^2/2)∑_m=0^M/21/m!(k^2σ^2/2)^m.

The summation is a truncated series expansion of exp(k^2σ^2/2) up to order M/2. This expression is the basis to prove that the HDAF filter approaches a true Dirac delta distribution in the limits of an infinitely broad filter or an infinitely large polynomial basis,

    lim_σ→ 0δ_M(k; σ) = lim_M→∞δ_M(k; σ) = 1   ∀ k∈ℝ,

thus identically preserving the function to reconstruct.

From equation (<ref>), one can also note that the Fourier expression of the HDAF is symmetric, bounded, and monotonically decreasing in k∈ℝ^+, with

    1 = δ_M(0; σ) ≥δ_M(k; σ) ≥lim_k→∞δ_M(k; σ) = 0,

therefore acting as a low-pass filter. Moreover, it has been shown that δ_M(k, σ) is an almost-ideal low-pass filter with transition frequency

    k^⋆ = √(M + 1) / σ,

and a transition region width that scales as O(M^-1/2σ^-1) <cit.>.



This behavior is depicted in Figure <ref> for varying M. There is a region of frequencies below k^⋆, the so-called DAF-plateau, such that δ_M(k; σ) ≈ 1. There follows a transition region centered around k^⋆ where the value of the filter smoothly vanishes, and then it indefinitely approximates to zero. As M is larger, the DAF plateau extends closer to k^⋆. Note that choosing σ from equation (<ref>) fixes the transition frequency k^⋆≈π/Δ x, which is the maximum frequency representable on a discrete uniform grid.

The HDAF reconstruction will be accurate for bandwidth-limited functions whose spectra lie within the DAF plateau, and any function with higher-frequency contributions will be smoothed. This suggests that this formalism is a good fit to use along with matrix product states since both techniques are especially suitable to represent bandwidth-limited functions <cit.>.




  §.§.§ Evaluation of the HDAF coefficients.



All the HDAF operators presented here are generated as a combination of displacements and coefficients

    K^(l)_τ = ∑_k=-2^n+1^2^n-1Δ x δ_M^(l)(kΔ x; σ, τ)(Σ^+)^k.

Moreover, the coefficients fulfill the general form

    Δ xδ_M^(l)(x; σ, τ)     = d_l∑_m=0^M/2h_m, l(x / √(2(σ^2 + iτ))),

with the definitions

    d_l   = (-1)^lΔ x/√(2(σ^2 + iτ))^l+1√(π),
    
      h_n, l(x)    = H_2n + l(x)exp(-x^2)c^n/n!, 
    
      c    = -1/4σ^2/(σ^2 + iτ).


From the properties of Hermite polynomials, it follows that h_n, l(x) obeys the double recurrence relation

    h_n+1, l(x)    = 2c/n + 1[ x h_n, l+1(x) - (2n + l + 1)h_n, l(x) ],
    
      h_n+1, l+1(x)    = 2x h_n+1, l(x) - 2c(2 + l/(n+1))h_n, l+1(x),

that makes the calculation of (<ref>) efficient and accurate, starting from the initial values h_0, l(x) = H_l(x)exp(-x^2) and h_0, l+1(x) = H_l+1(x)exp(-x^2).




  §.§.§ Effective summation bounds.



It was previously mentioned that equations (<ref>), (<ref>), (<ref>) and (<ref>) formally sum over all the grid points of x. In practice, however, the fast decay of the HDAF narrows the sum to a small subset of points around the origin. This subset can be further restricted since the filters are either symmetric or antisymmetric.

Only the highest power on the argument of δ_M^(l)(x; σ, τ) will contribute significantly to the value of the HDAF for x ≫ 0.
Let W be the smallest positive integer such that the coefficients contribute at most a predefined error tolerance ε,

    |Δ x δ_M^(l)(WΔ x; σ, τ)| ≤ε.


This integer can be estimated tightly by replacing the sum in (<ref>) with the highest power term on WΔ x from the polynomial on h_M/2, l(WΔ x / √(2(σ^2 + iτ))), but it leads to a transcendental equation for W. Instead, we approximate the sum by its last term as a whole and use the following upper bound for the Hermite polynomials with complex argument <cit.>,

    | H_n(z) | ≤√(2)^n√(n!)exp(√(2n)|z|),    z∈ℂ, n∈ℕ.


While this bound tends to overestimate the polynomial, the Gaussian envolvent will dominate fast enough for this not to become a significant problem. Then, we find W by setting

    |Δ xδ_M^(l)|    ≲Δ x√((M+l)!)/√(2π)√(|σ^2 + iτ|)^l+1(M/2)!|σ^2/2 (σ^2 + iτ)|^M/2
         ×exp( WΔ x√(M + l/|σ^2 + iτ|)) 
         ×exp(- W^2Δ x^2/2(σ^2 + τ^2/σ^2)) 
       = ε,

which reduces to the quadratic equation,

    W^2Δ x^2/2(σ^2 + τ^2/σ^2) - WΔ x√(M + l/|σ^2 + iτ|) + lnε/η = 0,

with


    η = Δ x√((M+l)!)/√(2π)√(|σ^2 + iτ|)^l+1(M/2)!|σ^2/2 (σ^2 + iτ)|^M/2.


The HDAF MPOs (<ref>) are obtained by summing over indices -W to W, where W is the closest integer from above to the solution of (<ref>). This sum contains 2W + 1 weighted displacement operators Σ̂^± k, where only W + 1 coefficients must be explicitly computed due to the symmetry δ_M^(l)(-x) = (-1)^lδ_M^(l)(x). Despite the number of summands, the resulting MPO will, in practice, be relatively simple, with a small bond dimension, for a reasonable choice of the HDAF parameters.




§ TIME EVOLUTION ALGORITHMS


Our goal is to solve the time-dependent Schrödinger equation (<ref>) with H=D(-∂_x^2)+V(x). The formal solution for this problem can be expressed as the repeated action of a possibly time-dependent unitary operator U(t) on an initial state ψ(x,t=0),

    ψ(x,t)=U(t)ψ(x,0)=e^-iĤtψ(x,0).


In the particular framework of problems we are interested in, the state ψ(x,t) will be encoded using MPS/QTT, and we will use MPOs and finite-precision MPS algebra tools to approximate the unitary operator U(t) for brief time periods. More explicitly, the studies below will either use an MPO structure to encode the Hamiltonian H or create an MPO that directly approximates U(t). In both cases, the representations based on QTT/MPS will provide us access to exponentially dense grids with 2^n points in space, which may prove advantageous regarding vector representations.

To address the time evolution problem, the PDE operators described in Section <ref>
require global evolution schemes independent of the locality of interactions. MPS algorithms present many suitable alternatives, such as the time-dependent variational principle (TDVP) <cit.>, and
Taylor, Padé, and Arnoldi approximations of the evolution operator <cit.>. This section presents a selection of time-evolution methods with an MPO-MPS implementation: explicit (Euler, Improved Euler, and fourth-order Runge-Kutta) and implicit (Crank-Nicolson) Runge-Kutta methods, restarted Arnoldi iteration, and the split-step method. All previous methods are suitable for a finite difference and HDAF approximation of the differential operator.
Additionally, HDAF's propagator approximation enables the use of the split-step method.



 §.§ Runge-Kutta methods

Runge-Kutta methods approximate the time evolution by a Taylor expansion of the state with a local error, i.e., one-step error, that scales algebraically with the expansion order m as O(Δ t^m+1). Let us describe some of the most representative variations.

1. Euler method.  The simplest order one method

    ψ_0    = ψ(x,t_0),  
    ψ_k+1   = ψ_k - iΔ t H ψ_k,    k=0,1,…,N-1.


2. Improved Euler or Heun method. This method improves on the Euler method with a second-order error given by

    ψ_k+1 = ψ_k - i Δ t/2[v_1 + H (ψ_k-iΔ t v_1)],

with v_1 = H ψ_k.

3. Fourth-order Runge-Kutta method. Finally, the well-known fourth-order scheme

    ψ_k+1   =ψ_k + iΔ t/6(v_1+2v_2+2v_3+v_4), 
    
      v_1    = - H ψ_k, 
    
      v_2    = - H(ψ_k+iΔ t/2v_1), 
    
      v_3    = - H(ψ_k+iΔ t/2v_2), 
    
      v_4    = - H(ψ_k+iΔ t v_3).

This is one of the most commonly used methods in solving PDEs due to its balance in accuracy, stability, and simplicity.

4. Crank-Nicolson method. Implicit methods can increase numerical stability. The Crank-Nicolson algorithm is a second-order implicit method based on the trapezoidal rule that combines the Euler method and its backward version evaluated on the k and k+1 iterations, respectively. Thus, the state at the k+1 iteration is approximated as

    (𝕀+iΔ t/2H)ψ_k+1=(𝕀-iΔ t/2H)ψ_k.

Matrix inversion methods may solve the system of equations in its matrix-vector implementation. Other approaches, such as conjugate gradient descent, can be extended for its implementation in an MPO-MPS framework.



 §.§ Restarted Arnoldi iteration

Another alternative is to use Krylov subspace methods; more concretely, the restarted Arnoldi iteration adapted to the time evolution problem. These methods rely on the Krylov basis 𝒦_L = lin{|ψ_k⟩, H|ψ_k⟩,…,H^L-1|ψ_k⟩} to construct an approximation of the evolution. 

The restarted Arnoldi iteration constructs a Krylov basis { v_i}_i=1,…,n_v of n_v elements and computes the matrices of the expectation value of the operator H and its norm—A and N, whose matrix elements are ⟨ v_i|H| v_j⟩ and ⟨ v_i| v_j⟩, respectively— to approximate the exact exponential evolution as

    ψ_k+1=e^-iΔ t N^-1Aψ_k.

The error in approximating the exponential function is limited by the number of Krylov vectors, which scales with O(Δ t^n_v). This means that even a low number of vectors—n_v=5,10—can provide a highly accurate approximation. As a result, the cost of matrix inversion and exponentiation decreases significantly compared to the exact application of the H operator since the dimensions of A and N can remain constant with the system size, thus avoiding exponential scaling.



 §.§ Split-step method

Split-step methods are based on an approximate decomposition of the Hamiltonian exponential as a product of exponentials that can be efficiently computed.

The first-order method relies on the Lie-Trotter product formula to approximate the evolution operator as

    U(Δ t) = e^-i Δ t D(-∂_x^2) e^-iΔ t V(x) + O(Δ t^2).

Higher order expansions, such as the Suzuki-Trotter formulas <cit.>, enable a better approximation by decreasing the error scaling with the time step. A common alternative is the second-order approximation

    e^-i Δ t (D(-∂_x^2)+V(x))   = e^-iΔ t V(x)/2 e^-i Δ t D(-∂_x^2)
         × e^-i Δ t V(x)/2 + O(Δ t^3).



This decomposition generates the Störmer-Verlet integration scheme, the most common and one of the simplest members of the class of symplectic integrators <cit.>. These integrators are designed to preserve the system's energy and are especially well-fitted for conservative long-time evolution problems. While most implementations rely on second-order accuracy in Δ t, higher-order schemes can be constructed at the expense of introducing a larger number of exponential operators <cit.>.

Usually, the split-step (<ref>) relies on two Fourier transforms to compute the action of the free propagator exp(-iΔ t D(-∂_x^2)), but the HDAF formalism permits the efficient alternative of applying the free-propagator approximation (<ref>) directly in the coordinate representation.

The potential propagator exp(-iΔ t/2V(x)), diagonal in the coordinate basis, is approximated using TT-cross interpolation as implemented in Ref. <cit.>.
While the propagator for the harmonic potential can be cast analytically as an MPO, other potentials such as Eq. (<ref>) do not have an exact representation and must be computed numerically.




 §.§ One-step study








All methods described before must be iteratively applied over small intervals Δt to compose the whole evolution of a quantum state ψ(x,t). A study and comparison of all methods over one such integration step is a good proxy for understanding the performance and accuracy of the algorithms over longer simulations.

The benchmark problem in this study will be the evolution of a quenched state under the Hamiltonian (<ref>). This analytically solvable problem can be used to estimate the errors in the wavefunction for different algorithms, grids, and time steps. Our simulation studies the expansion of a quantum state as we reduce the trapping frequency of the harmonic potential by a factor ω_H/ω_0=0.01. This leads to a 100-fold quantum state expansion, increasing its standard deviation from σ_0 to σ_max=100σ_0. The simulation domain in position space is designed to capture the wavepacket at its maximum expansion. This means we will represent functions over x ∈ [-L/2, L/2) with L = 16σ_max. Since the initial state is derived from a tightly confined potential, the initial wavefunction is narrowly concentrated around x=0, which sets a lower bound for the grid discretization and number of qubits to accurately represent the initial and final state. Note that loading a function that is mostly zero outside a narrow interval is a tricky task for the MPS, requiring us to pad the function with zeros to compensate for the possible sampling or representation errors of TT-Cross or Chebyshev methods on the tails of the exponential.

The accuracy and performance of the algorithms are gauged using three figures of merit: (i) the function norm-2 difference 

    ε =√(∑_i |ψ(x_i,Δ t)-ψ̃(x_i,Δ t)|^2Δ x),

measures the accuracy of the methods—ψ(x_i,Δ t) is the analytic solution for t=Δ t, and ψ̃(x_i,Δ t) is the one-step state approximated by the method—; and (ii) the run time and (iii) maximum bond dimension χ_max determine their cost. 



The one-step study must separately study the time and spatial discretizations to independently understand the performance of the numerical integration and the PDE representation methods. Regarding spatial discretization, section <ref> already demonstrated the exponential advantage of HDAF methods compared to finite differences in approximating the derivatives. With this idea in mind, the first study uses a grid with n=18 qubits (262144 points) to explore the different integration techniques. The finite difference approximation uses filter nine in Ref. <cit.> to enhance noise suppression, and we fix a value M=40 for the HDAF approximation. Both techniques use periodic boundary conditions.

Figure <ref> shows the error ε (<ref>) and the run time scaling with Δ t for finite difference and HDAF derivative approximation for all methods. The scaling has saturated to the theoretical one for each time evolution method using HDAF approximation (Figure <ref>(b)), achieving numerical precision only limited by the MPS truncation. However, the truncation error of the finite difference approximation limits the accuracy of the evolution in Figure <ref>(a), leading to an error above the plateau of the HDAF implementation. Appendix <ref> shows the exact numerical scaling of ε with Δ t for these numerical simulations.


As shown in Figures <ref>(c)-(d), the implementation cost of HDAF and finite-difference algorithms is similar. The explanation for this is that the cost is dominated by the bond dimensions of the MPS states, which are similar for both algorithms. Given this and the exponentially improved accuracy of HDAF operators, we will, from now on, focus on the HDAF spectral methods and abandon the finite difference approximations entirely.

All evolution algorithms involve a trade-off between accuracy—i.e., the approximation order—and the cost of the implementation, which Figures <ref>(b)-(d) accurately describe. In the HDAF implementation, all methods converge with the number of qubits in terms of the error ε. The methods show an algebraic dependence of the error ε with the step Δ t as predicted by the theory, except for small time steps, limited by the numerical accuracy of the MPS implementation. The Arnoldi n_v=10 method is the best-performing one in terms of accuracy since it enables the largest step sizes due to its higher order. The Runge-Kutta, Arnoldi n_v=5, and split-step methods reach similar accuracies for a considerable Δ t range. The smaller order of the Euler method and the conjugate gradient implementation in the Crank-Nicolson method limit the error obtained beyond the intrinsic MPS truncation error. 


To identify a method with a suitable cost-accuracy trade-off, it is essential to analyze the error ε in conjunction with the run time (Figure <ref>(d)). Despite the Arnoldi n_v=10's low error, its run time is two orders of magnitude slower than the split-step method. Consequently, even with smaller time steps Δ t, the split-step method demonstrates a better balance of cost and accuracy. Additionally, practical implementations often do not require such high accuracy, allowing for larger time steps. Thus, the split-step method with the HDAF approximation of the propagator is the optimal choice for studying the expansion problem.

Let us compare the performance of MPS methods to their vector implementation using split-step methods and the state-of-the-art fast Fourier transform (FFT). Figure <ref> analyzes the scaling of split-step MPS and vector-based implementations with the discretization size. The MPS method examines various values of the truncation tolerance for the SVD and the simplifications of the finite precision algebra, measured as the norm-2 difference ‖ψ-ϕ‖^2 between the original state |ψ⟩ with bond dimension χ_ψ and the one projected in the subspace of MPS with bond dimension χ_ϕ, MPS_χ_ϕ, such that χ_ϕ<χ_ψ.



Figure <ref>(a) demonstrates that the MPS tolerance dominates the split-step error. Achieving numerical precision comparable to the vector implementation requires tolerances of order O(10^-28) or smaller. This error scaling with the number of qubits demonstrates that a discretization with Δ x ≈ 10^-1 suffices for HDAF to converge well above the finite difference approach. Regarding the run time (Figure <ref>)(b)), the FFT is more efficient than the vector-based HDAF, and MPS perform asymptotically better than vector approaches, as they present an exponential scaling of time with the number of qubits due to the exponential increase in the number of points. The run time of the MPS method is similar regardless of the tolerance since all resulting states achieve a bond dimension of a similar order. Since tolerances below 10^28 overestimate the bond dimension needed, we choose this value for the SVD truncation in previous and future simulations while keeping a tolerance of the order of the numerical precision for MPS simplification. Finally, Figure <ref> shows that the MPS error ε scaling with the time step Δ t presents a similar behavior to the vector-based implementation for tolerances smaller than 10^-28. As the tolerance increases, it limits the accuracy for smaller time steps, for which the MPS truncation error dominates over the split-step truncation error. As the time step increases, the error associated with the method allows for larger tolerances.




§ QUANTUM QUENCH EVOLUTION


After discussing the MPS-based algorithms for solving time-dependent PDEs, this section demonstrates the utility of the methods in a physical application. The problem under study has been presented in Section <ref> and consists of the expansion of a quantum particle in a broad potential, which may be harmonic, as in Section<ref>, or present a double-well structure (Section <ref>). The study will focus on the algorithms that show the best balance between performance and accuracy, namely the split-step methods implemented with a fast-Fourier transform in the vector case and with HDAF in the MPS/QTT simulations.



 §.§ Harmonic expansion


The first set of simulations addresses a similar problem to Section <ref>, studying the expansion of a particle in a harmonic potential (<ref>) that is 100 times weaker than the trap that confines the initial wavepacket. The simulation results in a 100-fold expansion of the wavefunction, which is properly captured with a discretization using 2^20 = 1048576 points or n=20 qubits for x∈[L/2, L/2) with L=16σ_max. The expansion starts with the state of the original harmonic oscillator with frequency ω_0=1, which is a real-valued Gaussian function. The potential is instantaneously weakened, relaxing the trapping frequency to ω_H=0.01. This results in an acceleration of the wavepacket, an expansion that lasts until t_f=0.5π /ω_H, at which the Gaussian solution (<ref>) reaches its maximum width. Let us consider three values for the time steps Δ t = 0.01,0.1,1 and a final time t_f=158, so it fits all time steps. As in Section <ref>, x∈[L/2,L/2) with L=16σ_max.

Figures <ref>(a)-(b) present the scaling of the error ε (<ref>) with the time t of the evolution for the MPS HDAF and vector FFT implementations of the split-step method, respectively. Note how the errors follow algebraic laws with very similar coefficients (see Appendix <ref>) that we attribute to the Störmer-Verlet integration scheme used: The error growth in time is consistent with the linear accumulation expected from a symplectic integration algorithm on a periodic Hamiltonian system <cit.>, as opposed to the quadratic law expected from general non-symplectic integrators such as Runge-Kutta on the same scenario <cit.>. The saturation of the global error could be explained since symplectic integrators keep the energy errors bounded at all times <cit.>, also implying a bound on phase-space errors for conservative systems with periodic orbits. These properties make the split-step method very suitable for long-time simulations.


The run time scaling (Figures <ref>(c)-(d)) is close to linear for the FFT algorithms, as the problem size fixes the cost of each step, and the total run time is the sum of individual evolution steps. In contrast, the run time of the MPS simulation is slightly above linear, a fact that can be explained by the growth of the bond dimension as the wavefunction expands, leading to an increase in memory size and also in the cost of various MPS and MPO-MPS operations.





The MPS algorithms' run time depends on the state's bond dimension at each time. Figure <ref>(b) shows the maximum bond dimension χ_max for each time step. The numerical methods exhibit a similar behavior for the bond dimension of the solution in time, compared to the exact evolution, i.e., the bond dimension increases with the absolute value of the phase β(t) of the analytic solution (<ref>), with the lower bond dimensions associated to the initial and final states, which are real Gaussians. This behavior appears due to the chirping of the wavefunction, which is inherent to the physical setting and does not depend on the numerical method used, up to some precision allowance.
However, the bond dimension of the exact solution acted mainly as an upper bound to the bond dimension of the evolution we found, suggesting that the errors induced by our numerical methods deviated the solution not in a random direction but towards a more efficient MPS representation. This is consistent with the HDAF theory, where, before discretization, the approximations arise solely due to the attenuation of high frequencies. The HDAF operators seem to be well-suited to the MPS/QTT framework by relying on approximations that are a good fit for their formalism.


Figure <ref>(a) shows the pointwise error of the maximum width state for the particle expansion, |ψ(x,t)-ψ̃(x,t)|, where ψ̃(x,t) is the solution approximated by the numerical methods and  ψ(x,t) is the analytic solution (<ref>). We observe that both implementations of split-step have similar error shapes that differ on the extremes of the interval for the larger step size, possibly due to errors associated with the MPS representation.

Let us focus on a concrete case to study the evolution of the wavepacket. Figure <ref> depicts the evolution computed using the MPS HDAF split-step method for Δ t = 0.1. As predicted by the analytic solution (<ref>), the harmonic potential (Figure <ref>(a)) induces an expansion of the particle, which is depicted in Figure <ref>(b). 






 §.§ Double well potential
 

The calibration of simulation conditions from section <ref> allows us to make informed decisions on the split-step algorithm, spatial discretization, time steps, and truncation errors. Let us now use this information to discuss a more interesting problem: the expansion of a nanoparticle in an anharmonic double-well potential. This problem is equivalent to a “double-slit” experiment for the particle, which will ideally spread into a coherent superposition of both halves of the trapping potential, and it constitutes a similar application to other problems studied in levitodynamics <cit.>. For this particular simulation, the double-well potential is an open harmonic trap divided by a small Gaussian perturbation (<ref>), using u=1 and σ=1, and a trapping frequency that is once more 100 times smaller than the initial particle trap. The potential is depicted in Figure <ref>(a). As in the previous case, the harmonic term has a frequency ω_H<ω_0, which weakens the confinement and expands the particle. The Gaussian term is repulsive since u>0, separating the larger potential in two wells, with a barrier around x=0.

For this simulation, the final expansion time chosen is t_f=1000, which surpasses one period of the harmonic quantum quench solution T=π/ω_H. This larger number is designed to explore multiple cycles of expansion and contraction of the particle's wavefunction to reveal collapse and revival dynamics. Figure <ref> shows simulation results, both from the particle's dynamics perspective and the wavefunction's complexity. During the evolution, the wavefunction's density |ψ̃(x,t)|^2 (Figure <ref>(b)) evolves according to the interplay of both terms in the potential (<ref>). The Gaussian term induces a separation in the particle's probability density. In contrast, the harmonic term preserves the confinement. It determines the period of the evolution, where the time of the maximum harmonic expansion t_f=0.5π/ω_H≈ 157.1 coincides with the maximum spread of the particle. The new term modifies the behavior of the harmonic potential (<ref>), which is depicted in Figure <ref>(c). As expected, the harmonic potential induces cyclic expansion of the wavepacket, and the added Gaussian term leads to a barrier in the potential that divides the probability density into two localization peaks traveling in opposite directions.
The Gaussian barrier also modifies the behavior of the state's bond dimension. It is no longer cyclic, like in the harmonic case (Figure <ref>(d)), but instead seems to saturate, leading to a decrease in the linear scaling of the run time as the system evolves. 





§ CONCLUSIONS


This work has introduced an HDAF encoding of differential operators for PDEs in an MPS/QTT framework. This encoding has shown exponential accuracy and low resource scaling. 

The MPS HDAF encoding enables the design of quantum-inspired time evolution algorithms to solve time-dependent PDEs: explicit and implicit Runge-Kutta methods, restarted Arnoldi iteration, and a split-step method. In particular, the split-step method benefits from the approximate representation of the free propagator unitary operator, which standard finite difference schemes cannot efficiently approximate.

The time evolution methods combined with HDAF overcome their finite difference implementation in terms of accuracy while maintaining a similar cost. Additionally, the split-step method shows the best trade-off in accuracy and cost. The HDAF time evolution algorithms are also competitive with state-of-the-art vector representations, enabling exponentially efficient encodings of functions and moderate overheads in the simulation.

The expansion of a particle in a broad potential acts as a benchmark for the methods. This poses a challenging problem since it defies the MPS representation due to the appearance of a chirp and rapid oscillations in the phase. Despite this chirp, the MPO-MPS algorithm produces accurate results with moderate bond dimensions and adequate run time.

Note that while the MPS and FFT speeds are comparable, only the MPS algorithm can scale up these grid densities to more dimensions. We expect to further develop and optimize these routines for higher-dimensional problems, accelerating them through C/C++ backends and other low-level optimizations to improve our Python programs' run time prefactors.

The present implementation is based on the SElf-Explaining Matrix Product State (SeeMPS) library for Python [<https://github.com/juanjosegarciaripoll/seemps2>].




§ ACKNOWLEDGMENTS

The authors would like to thank Juan José Rodríguez-Aldavero for his help with implementing the TT-cross interpolation. This work has been supported by Spanish Projects No. PID2021-127968NB-I00 and No. PDC2022-133486-I00, funded by MCIN/AEI/10.13039/501100011033 and by the European Union “NextGenerationEU”/PRTR”1. PGM acknowledges the funding by “FSE invierte en tu futuro” through an FPU Grant FPU19/03590
and by MCIN/AEI/10.13039/501100011033. JJGR and PGM acknowledge support
from CSIC Interdisciplinary Thematic Platform (PTI) Quantum Technologies (PTIQTEP+).
JG was supported by the Chilean National Agency for Research and Development (ANID-Chile), program “Doctorado Nacional 2020”, scholarship No. 21202616.
The authors also gratefully
acknowledge the Scientific Computing Area (AIC), SGAI-CSIC, for their assistance
while using the DRAGO Supercomputer to perform the simulations.




§ AUTHOR CONTRIBUTION STATEMENT


JG developed the extension of HDAF to the MPS formalism, its theoretical and numerical study, and the code implementation. PGM introduced the HDAF operators in the time evolution schemes and conducted their numerical characterization in the one-step and long-time simulations. JJGR conceptualized this work and its research goals. JJGR and LT supervised the research, providing solutions to challenges encountered throughout the process. JG and PGM wrote the original draft, and JJGR reviewed and edited the manuscript.





§ ONE-STEP Ε SCALING WITH Δ T

Figures <ref>(a)-(b) show a fit for the error ε with Δ t for the methods in Section <ref>, for the finite difference and HDAF derivative approximation of the derivative, respectively. Tables <ref> and <ref> contain the concrete numerical data of the fit ε=CΔ t ^m. We use a piecewise linear fit and show the data for larger Δ t, which are not limited by the MPS accuracy.







§ HARMONIC QUANTUM QUENCH EVOLUTION SCALING


Figure <ref> shows the error ε scaling and run time of the harmonic quantum quench evolution with time. Tables <ref>-<ref> contain the concrete numerical data of the fits ε = Ct^m and T = Ct^m, where T is the run time.











