2023 A space–time DG method for the Schr¨odinger equation with variable potential ∗

We present a space–time ultra-weak discontinuous Galerkin discretization of the linear Schr¨odinger equation with variable potential. The proposed method is well-posed and quasi-optimal in mesh-dependent norms for very general discrete spaces. Optimal h -convergence error estimates are derived for the method when test and trial spaces are chosen either as piecewise polynomials, or as a novel quasi-Treﬀtz polynomial space. The latter allows for a substantial reduction of the number of degrees of freedom and admits piecewise-smooth potentials. Several numerical experiments validate the accuracy and advantages of the proposed method.


Introduction
In this work we are interested in the approximation of the solution to the time-dependent Schrödinger equation on a space-time cylinder Q T = Ω×I, where Ω ⊂ R d (d ∈ N) is an open, bounded polytopic domain with Lipschitz boundary ∂Ω, and I = (0, T ) for some final time T > 0: ψ(x, 0) = ψ 0 (x) on Ω. (1.1) Here i is the imaginary unit; ∂ nx (•) is the normal derivative-in-space operator; V : Q T → R is the potential energy function; ϑ ∈ L ∞ (Γ R × I) is a positive "impedance" function; the Dirichlet (g D ), Neumann (g N ), Robin (g R ) and initial condition (ψ 0 ) data are given functions; Γ D , Γ N , Γ R are a polytopic partition of ∂Ω.
The model problem (1.1) has a wide range of applications.In quantum physics [25], the solution ψ is a quantum-mechanical wave function determining the dynamics of one or multiple particles in a potential V .In electromagnetic wave propagation [24], it is called "paraxial wave equation" and ψ is a function associated with the field component in a two-dimensional electromagnetic problem where the energy propagates at small angles from a preferred direction.In such problems, the function V depends on the refractive index and the wave number.In underwater sound propagation [22], it is referred to as "parabolic equation" and ψ describes a time harmonic wave propagating primarily in one direction.In molecular dynamics [2], by neglecting the motion of the atomic nuclei, the Born-Oppenheimer approximation leads to a Schrödinger equation in the semi-classical regime.
Space-time Galerkin methods discretize all the variables in a time dependent PDE at once; this is in contrast with the method of lines, which combines a spatial discretization and a time-stepping scheme.Space-time methods can achieve high convergence rates in space and time, and provide discrete solutions that are available on the whole space-time domain.
The literature on space-time Galerkin methods for the Schrödinger equation is very scarce.In fact, the standard Petrov-Galerkin formulation for the Schrödinger equation, i.e., the analogous formulation to that proposed in [32] for the heat equation, is not inf-sup stable, see [14, Sect.2.2].In [20], Karakashian and Makridakis proposed a space-time method for the Schrödinger equation with nonlinear potential, combining a conforming Galerkin discretization in space and an upwind DG time-stepping.This method reduces to a Radau IIA Runge-Kutta time discretization in the case of constant potentials.Moreover, under some restrictions on the mesh that are necessary to preserve the accuracy of the method, it allows for changing the spatial mesh on each time-slab, but not for local time-stepping.A second version of the method, obtained by enforcing the transmission of information from the past through a projection, was proposed in [21].This version reduces to a Legendre Runge-Kutta time discretization in the case of constant potentials.Recently, some spacetime methods based on ultra-weak formulations of the Schrödinger equation have been designed.The well-posedness of such formulations requires weaker assumptions on the mesh.Demkowicz et al., in [8], the authors proposed a discontinuous Petrov-Galerkin (DPG) formulation for the linear Schrödinger equation.The method is a conforming discretization of an ultra-weak formulation of the Schrödinger equation in graph spaces.Well-posedness and quasi-optimality of the method follow directly from the inf-sup stability (in a graph norm) of the continuous Petrov-Galerkin formulation.In [14], Hain and Urban proposed a space-time ultra-weak variational formulation for the Schrödinger equation with optimal inf-sup constant.The formulation in [14] is closely related to the DPG method in [8], but differs in the choice of the test and trial spaces.While for the method in [8] one first fixes a trial space and then construct a suitable test space, the method in [14] requires the choice of a conforming test space and then the trial space is defined accordingly.We are not aware of publications proposing space-time DG methods for the Schrödinger equation other than [8,14,20,21], outlined in this paragraph, and the space-time Trefftz-DG method in [11,12], which motivated the present paper.
Trefftz methods are Galerkin discretizations with test and trial spaces spanned by local solutions of the considered PDE.Trefftz methods with lower-dimensional spaces than standard finite element spaces, but similar approximation properties, have been designed for many problems, e.g., Laplace and solid-mechanics problems [31]; the Helmholtz equation [16]; the time-harmonic [15], and time-dependent [10] Maxwell's equations; the acoustic wave equation in second-order [1] and first-order [27] form; the Schrödinger equation [11]; among others.Nonetheless, pure Trefftz methods are essentially limited to problems with piecewise-constant coefficients, as for PDEs with varying coefficients the design of "rich enough" finite-dimensional Trefftz spaces is in general not possible.A way to overcome this limitation is the use of quasi-Trefftz methods, which are based on spaces containing functions that are just approximate local solutions to the PDE.In essence, the earliest quasi-Trefftz spaces are the generalized plane waves used in [17] for the discretization of the Helmholtz equation with smoothly varying coefficients.More recently, a quasi-Trefftz DG method for the acoustic wave equation with piecewise-smooth material parameters was proposed in [19], where some polynomial quasi-Trefftz spaces were introduced.As an alternative idea, the embedded Trefftz DG method proposed in [23] does not require the local basis functions to be known in advance, as they are simply taken as a basis for the kernel of the local discrete operators in a standard DG formulation.This corresponds to a Galerkin projection of a DG formulation with a predetermined discrete space onto a Trefftz-type subspace.In practice, it requires the computation of singular or eigenvalue decompositions of the local matrices.
In [11], the authors proposed a space-time Trefftz-DG method for the Schrödinger equation with piecewise-constant potential, whose well-posedness and quasi-optimality in mesh-dependent norms were proven for general discrete Trefftz spaces.Optimal h-convergence estimates were shown for a Trefftz space consisting of complex-exponential wave functions.
In this work we propose a space-time DG method for the discretization of the Schrödinger equation with variable potentials, extending the formulation of [11] to more general problems and discrete spaces.The main advantages of the proposed method are the following: • The proposed ultra-weak DG variational formulation of (1.1) is well-posed, stable, and quasioptimal in any space dimension for an almost arbitrary choice of piecewise-defined discrete spaces and variable potentials.
• A priori error estimates in a mesh-dependent norm can be obtained by simply analyzing the approximation properties of the local spaces.
• The method naturally allows for non-matching space-like and time-like facets and all our theoretical results hold under standard assumptions on the space-time mesh, which make the method suitable for adaptive versions and local time-stepping.
• Building on [19], for elementwise smooth potentials, we design and analyze a quasi-Trefftz polynomial space with similar approximation properties of full polynomial spaces but with much smaller dimension, thus substantially reducing the total number of degrees of freedom required for a given accuracy.

Structure of the paper:
In Section 2 we introduce some notation on the space-time meshes to be used and the proposed ultra-weak DG variational formulation on abstract spaces.Section 3 is devoted to the analysis of well-posedness, stability and quasi-optimality of the method.In Sections 4.2 and 4.3 we prove optimal h-convergence estimates for the method when the test and trial spaces are taken as the space of piecewise polynomials or a novel quasi-Trefftz space, respectively.In Section 5 we present some numerical experiments that validate our theoretical results and illustrate the advantages of the proposed method.We end with some concluding remarks in Section 6.

Space-time mesh and DG notation
Let T h be a non-overlapping prismatic partition of Q T , i.e., each element K ∈ T h can be written as K = K x × K t for a d-dimensional polytope K x ⊂ Ω and a time interval K t ⊂ I.We use the notation . We call "mesh facet" any intersection F = ∂K 1 ∩ ∂K 2 or F = ∂K 1 ∩ ∂Q T , for K 1 , K 2 ∈ T h , that has positive d-dimensional measure and is contained in a d-dimensional hyperplane.We denote by n F = ( n x F , n t F ) ∈ R d+1 one of the two unit normal vectors orthogonal to F with n t F = 0 or n t F = 1.We assume that each internal mesh facet F is either a space-like facet if n x F = 0, or a time-like facet if n t F = 0.
We further denote the mesh skeleton and its parts as union of all the internal time-like facets, F space h := the union of all the internal space-like facets.
We employ the standard DG notation for the averages { {•} } and space • N and time • t jumps for piecewise complex scalar w and vector τ fields: where n x K ∈ R d and n t K ∈ R are the space and time components of the outward-pointing unit normal vectors on ∂K ∩ F time h and ∂K ∩ F space h , respectively.The superscripts "−" and "+" are used to denote the traces of a function on a space-like facet from the elements "before" (−) and "after" (+) the facet.
The space-time prismatic meshes described in this section may include hanging space-like and time-like facets, so the proposed method allows for local time-stepping and local space-time refinements.Tent-pitched meshes are popular in space-time methods for wave propagation problems; see e.g., [30] and [27,Eq. 3].However, such meshes do not lead to a semi-implicit discretization of the Schrödinger equation because the propagation speed of its solutions, which dictates the slope of space-like facets of the tents, is infinite.
We denote space-time broken function spaces as

Variational formulation of the DG method
For any finite-dimensional subspace V hp (T h ) of the broken Bochner-Sobolev space the proposed ultra-weak DG variational formulation for the Schrödinger equation (1.1) is: where for some mesh-dependent stabilization functions More conditions on these functions, in particular on their dependence on the local mesh size, will be specified in Section 4. The variational formulation (2.1) can be derived by integrating by parts twice in space and once in time in each element as in [11], and treating the Neumann and the Robin boundary terms similarly to [11,Rem. 3.7].However, as the current setting does not require the discrete space V hp (T h ) to satisfy the Trefftz property (Sψ |K = 0, ∀K ∈ T h ), there are an additional volume term that is needed to ensure consistency (the first integral over K in A (•; •)), and a local Galerkin-least squares correction term (the second integral over K in A (•; •)) that were not present in the previous method.Such additional terms vanish when V hp (T h ) is a discrete Trefftz space, thus recovering the formulation in [11].
Remark 1 (Implicit time-stepping through time-slabs).The variational problem (2.1) is a global problem involving all the degrees of freedom of the discrete solution for the whole space-time cylinder Q T .However, as upwind numerical fluxes are taken on the space-like facets, if the space-time prismatic mesh T h can be decomposed into time-slabs (i.e., if the mesh elements can be grouped in sets of the form Ω × [t n−1 , t n ] for a partition of the time interval of the form 0 = t 0 < t 1 < . . .< t N = T ), the global linear system stemming from (2.1) can be solved as a sequence of N smaller systems of the form for n = 2, . . ., N .This is comparable to an implicit time-stepping, and it naturally allows for local mesh refinement in different regions of the space-time cylinder Q T .Moreover, when T h is a tensor-product space-time mesh, the potential V does not vary in time, and the partition of the time interval is uniform, the matrices K n and R n are the same for every time-slab.
Remark 2 (Self-adjointness and volume penalty term).The well-posedness of the variational formulation (2.1) strongly relies on the L2 (K)-self-adjointness of the Schrödinger operator S(•) on each K ∈ T h (in the sense that K Sψ ϕ dV = K ψ Sϕ dV for all ψ ∈ V(T h ), ϕ ∈ C ∞ 0 (K), thanks to the fact that the only odd derivative in S is multiplied to the imaginary unit), which makes the local Galerkin-least squares correction term consistent.On the one hand, such term is essential in the proof of coercivity of the sesquilinear form A (• ; •) (see Proposition 1 below).On the other hand, numerical experiments suggest that it can be neglected without losing accuracy and stability, see Section 5.1.2below.This is also the case for the quasi-Trefftz DG method for the Helmholtz equation [18, §5.1.3]and for the wave equation [19, §5.1],where a similar correction term was used.Nonetheless, in the design of an ultra-weak DG discretization for a PDE with a non-self-adjoint differential operator L(•) (e.g., the heat operator L(•) = (∂ t − ∆ x ) (•)), the corresponding local least-squares correction term K∈T h K µLψ hp Lv hp dV would not control the consistency term K∈T h K ψ hp L * v hp dV arising from the integration by parts.
Remark 3 (Time-dependent potentials).The variational problem (2.1) allows for time-dependent potentials V .This is an important feature as, in such a case, the method of separation of variables cannot be used to reduce the time-dependent problem (1.1) to the time-independent Schrödinger equation.

Well-posedness, stability and quasi-optimality of the DG method
The theoretical results in this section are derived for any spatial dimension d, and are independent of the specific choice of the discrete space V hp (T h ).
Recalling that the volume penalty function µ, the stabilization functions α, β and the impedance function ϑ are positive, and that δ ∈ (0, 1  2 ), we define the following mesh-dependent norms on V(T h ): .
The sum of the L 2 (K)-type terms ensures that |||•||| DG + is a norm.That |||•||| DG is a norm on V(T h ) follows from the following reasoning (see also [11,Lemma 3.1]): if w ∈ V(T h ) and w DG = 0, then w is the unique variational solution to the Schrödinger equation (1.1) with homogeneous initial and boundary conditions.Moreover, by the energy conservation (if , for all t ∈ (0, T ]; therefore, w = 0.The DG norms in (3.1)-(3.2) are chosen in order to ensure the following properties of the sesquilinear form A (•; •) and the antilinear functional ℓ(•), from which the well-posedness and quasi-optimality of the method (2.1) follow.

Proposition 1 (Coercivity). For all w ∈ V(T h ) the following identity holds
Proof.The result follows from the following identities (see [11,Prop. 3.2] for more details): Proposition 2 (Continuity).The sesquilinear form A (•; •) and the antilinear functional ℓ(•) are continuous in the following sense: ) and F D h are controlled as in [11,Prop. 3.3].The remaining terms are bounded using Cauchy-Schwarz inequality and the inequality δ ≤ 1 − δ < 1.
Theorem 1 (Quasi-optimality).For any finite-dimensional subspace V hp (T h ) of V(T h ), there exists a unique solution ψ hp ∈ V hp (T h ) satisfying the variational formulation (2.1).Additionally, the following quasi-optimality bound holds: Moreover, if g D = 0 and g N = 0 (or Proof.Existence and uniqueness of the discrete solution ψ hp ∈ V hp (T h ) of the variational formulation (2.1), and the quasi-optimality bound (3.4) follow directly from Propositions 1-2, the consistency of the variational formulation (2.1) and Lax-Milgram theorem.The continuous dependence on the data (3.5)follows from Proposition 1, and the fact that if g D = 0 and g N = 0 (or Γ D = ∅ and Γ N = ∅), the term w DG + on the right-hand side of (3.3b) can be replaced by w DG .Theorem 1 implies that it is possible to obtain error estimates in the mesh-dependent norm ||| • ||| DG by studying the best approximation in V hp (T h ) of the exact solution in the ||| • ||| DG + norm.Moreover, according to Proposition 3 below, a priori error estimates can be deduced from the local approximation properties of the space V hp (T h ) only, as the ||| • ||| DG + norm can be bounded in terms of volume Sobolev seminorms and norms.The proof of error estimates in mesh-independent norms on the full computational domain for ultra-weak DG methods is a delicate issue; see e.g., [15,Lemma 1] and [27, §5.4] for related results concerning Trefftz methods for the Helmholtz and the wave equations, respectively.
So far, we have not imposed any restriction on the space-time mesh T h .Henceforth, in our analysis we assume: • Uniform star-shapedness: There exists 0 < ρ ≤ 1 2 such that, each element K ∈ T h is star-shaped with respect to the ball B := B ρhK (z K , s K ) centered at (z K , s K ) ∈ K and with radius ρh K .
• Local quasi-uniformity in space: there exists a number lqu(T The proof of Proposition 3 is a direct consequence of a collection of trace inequalities (see [3,Theorem 1.6.6]and [27, Lemma 2]), which in our space-time setting can be written for any element where D 2 x ϕ is the spatial Hessian of ϕ, and C tr ≥ 1 only depends on the star-shapedness parameter ρ.
2 ), and assume that V ∈ L ∞ (K), ∀K ∈ T h .For all ϕ ∈ V(T h ), the following bound holds , where The factor 3 2 C tr appearing in the bound of Proposition 3 is due to the integral terms with arguments The existence and uniqueness of the discrete solution for any choice of the discrete space V hp (T h ), as well as the quasi-optimality estimate (3.4), follow from the coercivity and continuity of the sesquilinear form A (•; •) on the continuous space V(T h ) in Propositions 1 and 2, together with the consistency of the method.Thus, optimal convergence rates can be proven for the full polynomial space as in Section 4.2, since this space provides a good enough approximation of any sufficiently smooth solution.On the other hand, the quasi-Trefftz space introduced in Section 4.3 would require some adjustments in order to approximate the solution of an inhomogeneous problem.
The proposed DG method is dissipative, but the energy loss can be quantified in terms of the local least-squares error, the initial condition error, the jumps of the solution on the mesh skeleton, and the error on F D h ∪ F N h due to the weak imposition of the boundary conditions.More precisely, for g D = 0, g N = 0 and F R h = ∅, the discrete solution to (2.1) satisfies where This follows from the definition of the ||| • ||| DG norm of the solution ψ hp , the coercivity of the sesquilinear form A (• ; •), the definition of the antilinear functional ℓ(•) and simple algebraic manipulations; see [11,Rem. 3.6].

Discrete spaces and error estimates
In this section we prove a priori h-convergence estimates on the ||| • ||| DG + norm of the error for some discrete polynomial spaces.In particular, for each element K ∈ T h , we consider two different polynomial spaces: the space P p (K) of polynomials of degree p on K, and a quasi-Trefftz subspace QT p (K) ⊂ P p (K) with much smaller dimension, i.e., dim(QT p (K)) ≪ dim(P p (K)) (see Proposition 5 below).A polynomial Trefftz space for the case of zero potential V has been studied in [12].We denote the local dimensions n d+1,p := dim(QT p (K)) and r d+1,p := dim(P p (K)) in dependence of the space dimension d of the problem and the polynomial degree p, but independent of the element K.For simplicity, we only describe the case where the same polynomial degree is chosen in every element; the general case can easily be studied.

Multi-index notation and preliminary results
We use the standard multi-index notation for partial derivatives and monomials, adapted to the space-time setting: for j = (j x , j t ) = (j x1 , . . ., j We also recall the definition and approximation properties of multivariate Taylor polynomials, which constitute the basis of our error analysis.On an open and bounded set Υ ⊂ R d+1 , the Taylor polynomial of order m ∈ N (and degree m − 1), centered at (z, s) ∈ Υ, of a function ϕ ∈ C m−1 (Υ) is defined as If ϕ ∈ C m (Υ) and the segment [(z, s), (x, t)] ⊂ Υ, the Lagrange's form of the Taylor remainder (see [4,Corollary 3.19]) is bounded as follows: where h Υ is the diameter of Υ.In particular, if Υ is star-shaped with respect to (z, s), then the following estimate is obtained which, together with the well-known identity (see [3, Prop.(4.1.17)]) The Bramble-Hilbert lemma provides an estimate for the error of the averaged Taylor polynomial, see [9] and [3,Thm. 4.3.8].
, be an open and bounded set with diameter h Υ , star-shaped with respect to the ball B := B ρhΥ (z, s) centered at (z, s) ∈ Υ and with radius ρh Υ , for some 0 < ρ ≤ 1 2 .If ϕ ∈ H m (Υ), the averaged Taylor polynomial of order m (and degree m − 1) defined as satisfies the following error bound for all s < m A sharp bound on C d,m,ρ > 0 is given in [9, p. 986] in dependence of d, s, m and ρ, and the second bound is proven in [27, Lemma 1].

Full polynomial space
In next theorem, we derive a priori error estimates for the DG formulation (2.1) for the space of elementwise polynomials Theorem 2. Let p ∈ N, fix δ as in Proposition 3 and assume that V ∈ L ∞ (Q T ).Let ψ ∈ V(T h )∩H p+1 (T h ) be the exact solution of (1.1) and ψ hp ∈ V hp (T h ) be the solution to the variational formulation (2.1) with V hp (T h ) given by (4.2).Set the volume penalty function and the stabilization functions as then the following estimate holds Moreover, if h Kx ≃ h Kt for all K ∈ T h , there exists a positive constant C independent of the element sizes h Kx , h Kt , but depending on the degree p, the L ∞ (Q T ) norm of V , the trace inequality constant C tr in (3.6), the local quasi-uniformity parameter lqu(T h ) and the star-shapedness parameter ρ such that Proof.The proof follows from the choice of the volume penalty function µ and the stabilization functions α, β, the quasi-optimality bound (3.4), Proposition 3, the inequality for all elements K ∈ T h , and the Bramble-Hilbert lemma 1.

Quasi-Trefftz spaces
We now introduce a polynomial quasi-Trefftz space.Let p ∈ N and assume that V ∈ C p−2 (K).For each K ∈ T h we define the following local polynomial quasi-Trefftz space: for some point (x K , t K ) in K.We consider the following global discrete space , then by the multi-index Leibniz product rule for multivariate functions we have where The next proposition is the key ingredient to prove optimal convergence rates in Theorem 3 for the DG method (2.1) when V hp (T h ) is chosen as the quasi-Trefftz polynomial space defined in (4.3).[ψ] at (x K , t K ) that appear in (4.5) are at most of total order |j|+ 2 ≤ p, so they coincide with the corresponding derivatives of ψ.Furthermore, since Sψ = 0, then ) be the exact solution of (1.1) and ψ hp ∈ V hp (T h ) be the solution to the variational formulation (2.1) with V hp (T h ) given by (4.4).Set the volume penalty function µ and the stabilization functions α, β as in Theorem 2.Then, the following estimate holds Moreover, if h Kx ≃ h Kt for all K ∈ T h , there exists a positive constant C independent of the mesh size h, but depending on the degree p, the L ∞ (Q T ) norm of V , the trace inequality constant C tr in (3.6), the local quasi-uniformity parameter lqu(T h ) and the measure of the space-time domain Proof.The proof follows from the choice of the volume penalty function µ and the stabilization functions α, β, the quasi-optimality bound (3.4), bound (3), the inequality The a priori error estimate in Theorem 3 requires stronger regularity assumptions on ψ than Theorem 2 (namely ψ ∈ C p+1 (T h ) instead of ψ ∈ H p+1 (T h )) due to the fact that QT p (K) is tailored to contain the Taylor polynomial T p+1 (xK ,tK ) [ψ], but in general it does not contain the averaged Taylor polynomial Q p+1 [ψ].
Remark 6 (Non-polynomial spaces).Optimal h-convergence estimates can also be derived for non-polynomial spaces, by requiring the local space V hp (K) to contain an element whose Taylor polynomial coincides with that of the exact solution.This is the approach in [11] for the Trefftz space of complex exponential wave functions for the Schrödinger equation with piecewise-constant potential.

Basis functions and dimension
So far, we have not specified the dimension and a basis for the space QT p (K), which is the aim of this section.
Recalling that r d,p = dim be bases of P p (R d ) and P p−1 (R d ), respectively.We define and the following n d+1,p elements of where g x K , • denotes the restriction of g : K , where x Any element q p ∈ QT p (K) can be expressed in the scaled monomial basis as for some complex coefficients {C j } |j|≤p .By the conditions D j Sq p (x K , t K ) = 0 for all |j| ≤ p − 2, in the definition of QT p (K), we have the following relations between the coefficients which can be rewritten as The conditions imposed in (4.6) on the restriction of b J to x 1 = x K fix the coefficients of their expansion for all j with j x1 ∈ {0, 1}.In Figures 1 and 2, we illustrate how the coefficients that are not immediately determined by the conditions in (4.6) (i.e., those for j x1 ≥ 2) are uniquely defined and can be computed for the (1 + 1)-and (2 + 1)-dimensional cases using the recurrence relation (4.7).Proposition 5.The set of functions {b J } n d+1,p J=1 defined in (4.6) are a basis for the space QT p (K).Therefore, Proof.We first observe that the set of polynomials {b J } n d+1,p J=1 is linearly independent due to their restrictions to x 1 = x (1) K .On the other hand, the relations (4.7), imply that q p is uniquely determined by its restriction q p (x (1) K , •) and the restriction of its derivative ∂ x1 q p (x (1) K , •).In addition, there exist some complex coefficients {λ s } n d+1,p s=1 such that case.The colored dots in the (j x , j t ) plane represent the coefficients C jx jt .Each shape connects three dots located at the points (j x , j t + 1), (j x , j t ) and (j x + 2, j t ): this shape represents one of the equations (4.7) which, given C jx(jt+1) and C jx jt , allows to compute C (jx+2)jt .If the 2p+1 values with j x ∈ {0, 1} (corresponding to the blue nodes in the shaded region) are given, then these relations uniquely determine all the other coefficients, which can be computed sequentially using the relations (4.7) by proceeding left to right in the diagram.In the figure p = 7, the number of nodes is r 2,p = 36, the number of nodes in the shaded region is n 2,p = 15, the number of relations is r 2,p − n 2,p = 21.
whence q p = n d+1,p s=1 λ s b s , which completes the proof.

Remark 7 (Quasi-Trefftz basis construction: difference between Schrödinger and wave equations).
The definition of the basis functions b J in (4.6) can be modified by fixing the restriction of b J and its partial derivative ∂ x ℓ b J to x ℓ = x (ℓ) K for any 1 ≤ ℓ ≤ d.However, it is not possible to assign the values for a given time t = t K , as the order of the time derivative appearing in the Schrödinger equation is lower than the order of the space derivatives.How this affects the basis construction is visible from Figure 1: the coefficients (the colored dots) can be computed sequentially when all the other coefficients of a relation (the Y-shaped stencil) are known, so it is possible to reach all dots moving left to right, but not moving bottom to top.Imposing the values at a given time is possible for the wave equation, as it is done in [19, §4.4], precisely because in that case time and space derivatives have the same order.
Remark 8 (Constant-potential case).The space QT p (K) does not reduce to a Trefftz space for the case of constant potential V .Nonetheless, the pure Trefftz space T p (K) defined as does not possess strong enough approximation properties to guarantee optimal h-convergence.In particular, it does not contain the Taylor polynomial of all local solutions to the Schrödinger equation; for d = 1, p = 1 and V = 0, T p (K) = span {1, x}; however, ψ(x, t) = exp x + i 2 t satisfies Sψ = 0, and T p+1 (0,0 Remark 9 (Trefftz dimension).As seen in Proposition 5, the quasi-Trefftz polynomial space has considerably lower dimension than the full polynomial space of the same degree.This "dimension Figure 2: A representation of the relations defining the coefficients of b J for the (2+1)-dimensional case.The colored dots in position j = (j x , j y , j t ), |j| ≤ p, correspond to the coefficients C jx jy jt (here p = 5 and r p = 56).Each white circle is connected by the segments to four nodes and represents one of the equations in (4.7): given C jx jy jt , C jx jy(jt+1) and C jx(jy +2)jt , it allows to compute C (jx+2)jy jt (the leftmost of the four nodes connected to a given white circle) using (4.7).
reduction" is common to all Trefftz and quasi-Trefftz schemes.In particular, the dimension n d+1,p of QT p (K) is equal to the dimension of the space of harmonic polynomials of degree ≤ p in R d+1 , the Trefftz space of complex exponential wave functions for the Schrödinger equation with piecewiseconstant potential in [11], the Trefftz and quasi-Trefftz polynomial space for the wave equation in [27, Eq. ( 42)-( 43)] and [19].

Numerical experiments
In this section we validate the theoretical results regarding the h-convergence of the proposed method, and numerically assess some additional features such as p-convergence and conditioning.Although we do not report the results here, optimal convegence rates of order O h p+1 are observed for the error in the L 2 (Q T )-norm.We list some aspects regarding our numerical experiments • We use Cartesian-product space-time meshes with uniform partitions along each direction, which are a particular case of the situation described in Remark 1.
• We choose (x K , t K ) in the definition of the quasi-Trefftz space QT p (K) in (4.3) as the center of the element K.
• In all the experiments we consider Dirichlet boundary conditions.
• The linear systems are solved using Matlab's backslash command.
• The quasi-Trefftz basis functions {b J } n d+1,p J=1 are constructed by choosing m J and m J in (4.6) as scaled monomials and by computing the remaining coefficients C j with the relations (4.7).
• In the h-convergence plots, the numbers in the yellow rectangles are the empirical algebraic convergence rates for the quasi-Trefftz version (continuous lines).The dashed lines correspond to the errors obtained for the full polynomial space.

(1 + 1)-dimensional test cases
We first focus on the (1 + 1)-dimensional case, for which families of explicit solutions are available for some well-known potentials V .

h-convergence
In order to validate the error estimates in Theorems 2 and 3, we consider a series of problems with different potentials V .No significant difference in terms of accuracy between the quasi-Trefftz and the full polynomial versions of the method with the same polynomial degree p (corresponding to different numbers of DOFs n d+1,p and r d+1,p , respectively) is observed in all the experiments.
In Figure 4, we present the errors obtained for ω = 10, n = 2 and a sequence of Cartesian meshes with uniform partitions and h x = h t = 0.05 × 2 −i , i = 0, . . . 4. Rates of convergence of order O (h p ) in the DG norm are observed, as predicted by the error estimate in Theorem 3. A convergence of at least order O h p+1 is observed for the L 2 -error at the final time, which is faster (by a factor h) than the order that can be deduced from the estimates in Theorems 2 and 3. We have also included the plots for the error decay with respect to the total number of degrees of freedom, where the same h-convergence rates are observed for both versions of the method (see also the p-convergence plot in Figure 9a for a clearer understanding of the dependence of the error on p).
(a) Energy error evolution  Due to the fast decay of the exact solution close to the boundary (see Figure 8 (panel a), the energy is expected to be preserved.In Figure 3, we show the evolution of the energy error, and the convergence of the energy loss E loss to zero for the quasi-Trefftz version.In the latter, rates of order O h 2p are observed, which follows from Remark 5 and the error estimates in Theorems 2 and 3. (a) DG error     ) quantum harmonic oscillator problem with potential V (x) = 50x 2 and exact solution ψ 2 in (5.1).Convergence with respect to the mesh size h (top panels) and the total number of degrees of freedom (bottom panels).
Reflectionless potential (V (x) = −a 2 sech 2 (ax)) This potential was studied in [5] as an example of a reflectionless potential.On the space-time domain Q T = (−5, 5) × (0, 1), we consider the Schrödinger equation with exact solution (see [13,Problem 2.48 In Figure 5, we show the errors obtained for a sequence of meshes with h x = 2h t = 0.2 × 2 −i , i = 0, . . ., 4, and a = 1.As in the previous experiment, rates of convergence of order O (h p ) and O h p+1 are observed in the DG norm and the L 2 norm at the final time, respectively.The real part of the exact solution is depicted in Figure 8 (panel b).
Morse potential (V (x) = D(1 − e −αx ) 2 ) This potential was introduced by Morse in [28] to obtain a quantum-mechanical energy level spectrum of a vibrating, non-rotating diatomic molecule.There, the following family of solutions was presented (see also [6]) where ⌊•⌋ is the floor function, n = 0, . . ., ⌊λ − 1/2⌋, L denote the general associated Laguerre polynomials as defined in [29,Table 18.3.1]and In Figure 6, we show the errors obtained for the Morse potential problem with D = 8, α = 4 and exact solution ψ 1,1 on the space-time domain Q T = (−0.5,1.5) × (0, 1) for a sequence of meshes with h x = h t = 0.1 × 2 −i , i = 0, . . ., 4. The observed rates of convergence are in agreement with those obtained in the previous experiments.The real part of the exact solution is depicted in Figure 8 (panel c).

Square-well potential
We now consider a problem taken from [11], whose exact solution is not globally smooth.On the space-time domain 2) × (0, 1), we consider the Schrödinger equation with homogeneous Dirichlet boundary conditions and the following square-well potential for some fixed V * > 0. The initial condition is taken as an eigenfunction (bound state) of where k * is a real root of the function f (k ).The solution of the corresponding initial boundary value problem (1.1) is ψ(x, t) = ψ 0 (x) exp(−ik 2 t) and belongs to the space H p+1 (T h ) ∩ C ∞ I; C 1 (Ω) \C ∞ I; C 2 (Ω) for all p ∈ N, provided that T h is aligned with the discontinuities of the potential V ; therefore, Theorems 2 and 3 apply.Among the finite set of values k * for a given V * , in this experiment we take the largest one, corresponding to faster oscillations in space and time.
In Figure 7, we show the errors obtained for V * = 20 (k * ≈ 3.73188) and a sequence of meshes with h t = √ 2h x = 0.1 × 2 −i , i = 0, . . ., 4. Optimal convergence in both norms is observed for the errors of the quasi-Trefftz version of the method.

Effect of stabilization and volume penalty terms
In this experiment we are interested in the effect of neglecting some of the terms in the variational formulation (2.1).To do so, we consider the (1+1)-dimensional quantum harmonic oscillator problem with exact solution (5.1).In Tables 1-2 (quasi-Trefftz space) and 3-4 (full polynomial space) we present the errors in the DG-norm obtained for the same sequence of meshes and approximation degrees as in the previous section, for different combinations of the stabilization terms α, β and the volume penalty parameter µ.Although the proof of well-posedness of the method (2.1) relies on the assumption that α, β and µ are strictly positive, in our numerical experiments, the matrices of the arising linear systems are non-singular and optimal convergence rates are observed even when all these parameters are set to zero.Moreover, the errors obtained when α = 0 or β = 0 are smaller as some terms in the definition (3.1) of • DG vanish, while the presence of µ seems to have just a mild effect in the results.Not shown here, similar effects were observed for the error in the L 2 (F T h )-norm.

p-Convergence
We now study numerically the p-convergence of the method, i.e., for a fixed space-time mesh T h , we study the errors when increasing the polynomial degree p.We consider the (1 + 1)-dimensional problems above with the same parameters and the coarsest meshes for each case.In Figure 9, we compare the errors obtained for the method with the two choices for the discrete space V hp (T h ) analyzed in the previous sections: the full polynomial space (4.2) and the quasi-Trefftz polynomial  [17,11,1,30] but no proof is available yet (differently from the stationary case, [16, §3]).In general, for a (d+1)-dimensional problem, we expect exponential convergence of order O(e −b d √ NDoF s ) and O(e −c (d+1) √ NDoF s ) for the quasi-Trefftz and full-polynomial versions, respectively.

Conditioning
We now assess the conditioning of the stiffness matrix.In Figure 10 we compare the 2-condition number κ 2 (•) for the stiffness matrix K n defined in Remark 1, for the free particle problem V = 0 on the space-time domain Q T = (0, 1) × (0, 1).We consider the proposed polynomial quasi-Trefftz space in (4.4), the full-polynomial space in (4.2) and the pure-Trefftz space of complex exponential wave functions T p (T h ) proposed in [11].A basis {φ ℓ } 2p+1 ℓ=1 ⊂ T p (T h ) was defined in [11] as (5.5) We consider two choices for the parameters κ ℓ : the arbitrary choice used in [11] κ ℓ = −p, . . ., p, and the choice κ ℓ = 2πℓ/h x which makes the basis orthogonal in each element.The conditioning number κ 2 (K) for the quasi-Trefftz space, the full polynomial space, and the Trefftz space with orthogonal basis asymptotically grows as O h −1 for all p ∈ N, while for the Trefftz space with a non-orthogonal basis, asymptotically grows as O h −(2p+1) .Unfortunately, with higher dimensions and non-Cartesian elements, choosing the parameters and directions defining the basis functions {φ ℓ } so as to obtain an orthogonal basis is more challenging.

(2 + 1)-dimensional test cases
We now present some numerical test for space dimension d = 2.We recall that we use Cartesian space-time meshes with uniform partitions along each direction.

h-convergence
Singular time-independent potential (V (x, y) = 1 − 1/x 2 − 1/y 2 ) We consider the (2 + 1)dimensional problem on Q T = (0, 1) 2 × (0, 1) with exact solution (see [33]) ψ(x, y, t) = x 2 y 2 e it . (5.6) In Figure 11, we show the errors obtained for a sequence of meshes with h x = h y = h t = 0.1, 0.0667, 0.05, 0.04 and different degrees of approximation p.As in the numerical results for the (1 + 1)dimensional problems, we obtain rates of convergence of order O (h p ) in the DG norm, and O h p+1 in the L 2 norm at the final time.
We now consider a manufactured problem with a time-dependent potential (see [7]).On the space-time domain Q T = (0, 1) 2 × (0, 1) the exact solution is ψ(x, y, t) = ie i(t−1/2) 4 sech(x)sech(y). (5.7) In Figure 12 we show the errors obtained for the sequence of meshes from the previous experiment, and optimal convergence is observed in both norms.

p-convergence
In Figure 13 we show the results obtained for the p-version of the method applied to the (2 + 1)dimensional problems above, on the coarsest mesh.As expected, for the (2 + 1)-dimensional case, the error of the quasi-Trefftz version decays root-exponentially as O e −b √ N dof s .

Concluding remarks
We have introduced a space-time ultra-weak discontinuous Galerkin discretization for the linear Schrödinger equation with variable potential.The DG method is well-posed and quasi-optimal in mesh-dependent norms for any space dimension d ∈ N, and for very general prismatic meshes and discrete spaces.We proved optimal h-convergence of order O (h p ), in such a mesh-dependent norm, for two choices of the discrete spaces: the space of piecewise polynomials, and a novel quasi-Trefftz polynomial space with much smaller dimension.When the space-time mesh has a

Figure 1 :
Figure 1: A representation of the relations defining the coefficients of b J for the (1+1)-dimensional Energy loss at T = 1

Figure 3 :
Figure 3: Time-evolution of the energy error for the quantum harmonic oscillator problem with potential V (x) = 50x 2 and exact solution ψ 2 in (5.1).

Figure 4 :
Figure4: h-convergence for the (1 + 1) quantum harmonic oscillator problem with potential V (x) = 50x 2 and exact solution ψ 2 in (5.1).Convergence with respect to the mesh size h (top panels) and the total number of degrees of freedom (bottom panels).

Figure 8 :
Figure 8: Real part of the exact solutions for the (1 + 1)-dimensional problems.

2
and exact solution ψ 2 in (5.1) for different combinations of the stabilization parameters α, β and volume penalty parameter µ = 0.

Figure 10 :
Figure 10: Conditioning of the stiffness matrix for the DG method with different discrete spaces.
L 2 error at T = 1
{w} }| 2 on F time h in the definition (3.1) of the ||| • ||| DG norm.The Remark 4 (Inhomogeneous Schrödinger equation).The space-time ultra-weak DG variational formulation in (2.1) can be easily extended to approximate the solution to inhomogeneous Schrödingertype problems with a sufficiently smooth term f : Q T → C at the right-hand side of the first equation in (1.1); see [26, Ch. 3, § 10] for the well-posedness of such problems.In order to preserve the consistency of the method, it is necessary to add the following term to the antilinear functional ℓ(•):

Table 1 :
h Fx DG error Rate DG error Rate DG error Rate DG error Rate h-convergence for the quasi-Trefftz version applied to the quantum harmonic oscillator problem with potential V (x) = 50x 2 and exact solution ψ 2 in (5.1) for different combinations of the stabilization parameters α, β and volume penalty parameter µ = 0.

Table 2 :
h Fx DG error Rate DG error Rate DG error Rate DG error Rate h-convergence for the quasi-Trefftz version applied to the quantum harmonic oscillator problem with potential V (x) = 50x 2 and exact solution ψ 2 in (5.1) for different combinations of the stabilization parameters α, β and volume penalty parameter µ = 0.