NONTENSORIAL GENERALISED HERMITE SPECTRAL METHODS FOR PDES WITH FRACTIONAL LAPLACIAN AND SCHRÖDINGER OPERATORS

Abstract. In this paper, we introduce two families of nontensorial generalised Hermite polynomials/functions (GHPs/GHFs) in arbitrary dimensions, and develop efficient and accurate spectral methods for solving PDEs with integral fractional Laplacian (IFL) and/or Schrödinger operators in R. As a generalisation of the G. Szegö’s family in 1D (1939), the first family of multivariate GHPs (resp.GHFs) are orthogonal with respect to the weight function |x|2μe ́|x| 2 (resp. |x|) in R. We further construct the adjoint generalised Hermite functions (A-GHFs), which have an interwoven connection with the corresponding GHFs through the Fourier transform, and are orthogonal with respect to the inner product ru, vsHspRdq “ pp ́Δqs{2u, p ́ΔqvqRd associated with the IFL of order s ą 0. As an immediate consequence, the spectral-Galerkin method using A-GHFs as basis functions leads to a diagonal stiffness matrix for the IFL (which is known to be notoriously difficult and expensive to discretise). The new basis also finds remarkably efficient in solving PDEs with the fractional Schrödinger operator: p ́Δq ` |x| with s P p0, 1s and μ ą ́1{2 in R. We construct the second family of multivariate nontensorial Müntz-type GHFs, which are orthogonal with respect to an inner product associated with the underlying Schrödinger operator, and are tailored to the singularity of the solution at the origin. We demonstrate that the Müntz-type GHF spectral method leads to sparse matrices and spectrally accurate solution to some Schrödinger eigenvalue problems.

polynomials (GHPs)), through a second-order differential equation in an exercise problem. The GHPs are orthogonal with respect to the weight function | | 2 e´2 in R. Chihara was among the first who systematically studied the properties of the GHPs, and the associated generalised Hermite functions (GHFs): p p q p q :" e´2 {2 p q p q (orthogonal with respect to the weight function | | 2 in R) in his PhD thesis [10] entitled as "Generalised Hermite Polynomials" (1955). Later, some standard properties were collected in his book [11] (1978). Whereas the usual Hermite polynomials/functions are well-studied especially in spectral approximations, there have been very limited works on this generalised family (see, e.g., [30,[35][36][37] for the properties or further generalisations). In particular, the generalised Hermite spectral method in terms of algorithms, analysis and applications are still under-explored, which is indeed a topic worthy of investigation.
The main purpose of this paper is to construct two families of nontensorial GHPs/GHFs in arbitrary dimensions, and explore their applications in solving PDEs with the IFL and/or Schrödinger operators.
Firstly, we construct the -dimensional GHPs t , ,ℓ p qu (cf. (2.12)) and GHFs t p , ,ℓ p qu (cf. (2.13)), which are orthogonal with respect to the weight functions | | 2 e´| | 2 and | | 2 in R with ą´1 2 , respectively. In one dimension, they reduce to Szegö's GHPs/GHFs (up to a constant multiple). More importantly, we further introduce a family of adjoint generalised Hermite functions (A-GHFs) t q , ,ℓ p qu (cf. (2.25)) and derive some appealing properties that are essential for developing fast and accurate spectral algorithms. We show that the adjoint pair is closely interwoven through the Fourier transform F r p , ,ℓ sp q " i`2 q , ,ℓ p q, F r q , ,ℓ sp q " p´iq`2 p , ,ℓ p q. (1.1) Notably, by construction, the A-GHFs are orthogonal with respect to the inner product that induces the so-called Gagliardo semi-norm of the fractional Sobolev space pR q for P p0, 1s, that is, r q , ,ℓ , q , , s pR q "`p´∆q 2 q , ,ℓ , p´∆q 2 q , ,˘R " ℓ , (1.2) where p´∆q is the integral fractional Laplacian operator (cf. (2.28) and (2.29)). An immediate implication is that the use of A-GHFs as basis functions in the spectral-Galerkin approximation of the IFL leads to a diagonal stiffness matrix. In contrast, it has been a nightmare for computing this matrix in a usual tensorial Hermite spectral method when " 3 (cf. [29]). Moreover, this new basis offers efficient spectral algorithm for solving PDEs with the fractional Schrödinger operator: p´∆q`p q with p q " | | 2 or more general p q " | | 2 p q (where is smooth function of Schwartz class) with P p0, 1s and ą´1{2. In light of the orthogonality (1.2), the stiffness matrix under the Galerkin framework using the basis t q , ,ℓ u becomes diagonal, while the singular potential | | 2 can be treated as the weight function, since q , ,ℓ can be represented as a linear combination of t p , , u (cf. (2.20) and (2.25)). Indeed, the A-GHFs can provide a viable tool for the solutions of fractional Schrödinger problems (see, e.g., [5,6,22,47]).
It is noteworthy that the 3D GHPs (with " 0 and an appropriate scaling) reduce to the Burnett polynomials [8] (1936), which are mutually orthogonal with respect to the Maxwellian ℳp q " p2 q´3 {2 e´| | 2 {2 , and have proven to be a useful basis in solving kinetic equations (cf. [9,20] and the references therein). It is important to point out that the GHFs with " 0 are eigenfunctions of the harmonic oscillator (cf. (2.23)):´∆`| In fact, such an attractive property (in 2D) has been explored in [4] for computing the ground states and dynamics of the Bose-Einstein condensation. It is of fundamental and practical interest to search for the explicit eigen-functions of the Schrödinger operator with a more general potential or some variance of the operator in (1.3), which is the second purpose of this paper. The main finding (cf. Thm. 4.2) is that for ą maxp1´{2, 0q, there exists a family of Müntz-type GHFs t p ℋ , ,ℓ u (cf. (4.3)) satisfying´∆`2 Optimal basis for the Schrödinger operator: for , in (4. 16) In particular, for " 1{2, we find´∆´``p´1 With a proper scaling, this gives the eigen-pairs of the Schrödinger operator with the Coulomb potential: where is a nonzero constant (cf. Cor. 4.3). By construction, this new family t p ℋ , ,ℓ u in the radial direction turns out to be some special Müntz functions, so it is dubbed as Müntz-type for distinction. We remark that a Müntz polynomial ř "0 is generated by a Müntz sequence: 0 ă 1 ă 2 ă¨¨¨ă (cf. [31] (1914)), and the set of Müntz polynomials with 0 " 0, and real coefficients t u are dense in the space of continuous functions if and only if ř 8 "0´1 " 8 (cf. [7] and the references therein). Such a tool finds very effective in approximating singular solutions (see, e.g., [19,39]). We shall demonstrate in Section 4 that the Müntz-type GHF spectral-Galerkin approach is the method of choice of the Schrödinger eigenvalue problems with the fractional power potential in terms of both the efficiency and accuracy. In particular, the spectral accuracy can be achieved by using such basis function to match the singular behaviours of the eigenfunctions.
In Table 1, we provide a roadmap of two types of generalisations and some of their properties that are essential for developing efficient spectral algorithms for PDEs with integral fractional Laplacian in Section 2 and the Schrödinger eigenvalue problems in Section 4. In Subsection 2.4, we highlight the differences and connections with the relevant existing generalisations, and further testify that most of our constructions herein are new.

Multivariate nontensorial generalised Hermite polynomials/functions
In this section, we first make necessary preparations by introducing some notation and properties of the spherical harmonic functions. We then define the GHPs and GHFs upon the generalised Laguerre polynomials and spherical harmonics, and further construct a family of adjoint GHFs. We present various appealing properties of these basis functions, which are essential for the efficient spectral algorithms to be developed in the forthcoming section. We conclude this section by elaborating on their differences and connections with the most relevant Hermite-related polynomials/functions in literature.

Preliminaries
Let R " p´8, 8q, N " t1, 2,¨¨¨u, and N 0 :" t0u Y N. For P N, we denote by R the -dimensional Euclidean space with the inner product and norm x , y :" ř "1 , and " | | :" a x , y, respectively, for any , P R . Denote the unit vector along any nonzero vector byˆ" {| |. Recall the -dimensional spherical coordinates 1 " cos 1 ; 2 " sin 1 cos 2 ;¨¨¨¨¨¨;´1 " sin 1¨¨¨s in´2 cos´1; " sin 1¨¨¨s in´2 sin´1, 1 ,¨¨¨,´2 P r0, s,´1 P r0, 2 s, with the spherical volume element In spherical coordinates, the -dimensional Laplacian takes the form where ∆ S´1 is the Laplace-Beltrami operator on the unit sphere S´1 :" t P R : | | " 1u. Define the inner product of 2 pS´1q as We next introduce the -dimensional spherical harmonics as in [12]. Let be the space of all realdimensional homogeneous polynomials of degree as follows As an important subspace of , the space of all real harmonic polynomials of degree is defined as It is known that the dimensionality dimp q "ˆ`´1˙, dimpℋ q "ˆ`´1˙´ˆ`´3 2˙: where it is understood that for " 0, 1, the value of the second binomial coefficient is zero (cf. [12] (1.1.5)). In fact, for " 1, all harmonic polynomials are spanned by t1, u.
The -dimensional spherical harmonics are the restrictions of harmonic polynomials in ℋ to S´1, denoted by ℋˇˇS´1. It is important to remark the correspondence between a harmonic polynomial and the related spherical harmonic function (cf. [12], Ch. 1): for any p q P ℋ , with pˆq P ℋˇˇS´1 . It is noteworthy that p q is a homogeneous polynomial in R , while pˆq is a nonpolynomial function on the unit sphere. For P N 0 , let t ℓ : 1 ď ℓ ď u be the real (orthogonal) spherical harmonic basis of ℋ | S´1 , and note that the spherical harmonics of different degree are mutually orthogonal (cf. [12], Thm. 1.1.2), i.e., ℋ | S´1 K ℋ | S´1 for " . Thus, we have where we introduce two-related the index sets Υ 8 " pℓ, q : 1 ď ℓ ď , 0 ď ă 8, ℓ, P N 0 ( , The spherical harmonic basis functions are eigenfunctions of the Laplace-Beltrami problem ∆ S´1 ℓ pˆq "´p`´2q ℓ pˆq. (2.10) The second building block of the GHPs/GHFs is the generalized Laguerre polynomials, denoted by p q p q for P p0, 8q and ą´1. They are orthogonal with respect to the weight function e´, that is, where we refer to [43], ( [25], Ch. 4) and ( [38], Ch. 7) for more properties.

Generalized Hermite polynomials/functions in R
Definition 2.1. For ą´1 2 , P N 0 and pℓ, q P Υ 8 , we define the -dimensional generalised Hermite polynomials as , ,ℓ p q :" , ,ℓ p ; q " p`´2 2`q p| | 2 q ℓ p q " p`´2 2`q p 2 q ℓ pˆq, "ˆ, (2.12) and the -dimensional generalised Hermite functions as (2.13) Remark 2.2. As we shall see later, the one-dimensional GHPs (up to a constant multiple) coincide with the one-dimensional generalisation first introduced in Szegö ( [43], p. 371) (1939), after which we name the above new families. Indeed, they include several special types of multivariate Hermite polynomials/functions with many applications in both theory and numerics. For example, the three-dimensional GHPs with " 0 and a proper scaling lead to the Burnett polynomials [8] (1936) which have rich applications in kinetic theory (see [9] and the references therein). The notion of constructing special Laguerre-Fourier basis functions (relevant to the two-dimensional GHPs with " 0) for computing the ground states and dynamics of Bose-Einstein condensation [34] was found effective in e.g., [4].
Before we consider the applications of GHPs and GHFs, we first present some attractive properties. By construction, they enjoy the following important orthogonality. Theorem 2.3. For ą´1 2 , , P N 0 and pℓ, q, p , q P Υ 8 , the GHPs are mutually orthogonal with respect to the weight function | | 2 e´| | 2 , namely, 14) and the GHFs are orthonormal, viz., Proof. The orthogonality (2.15) is a direct consequence of (2.13) and (2.14), so we only need to show (2.14).
In view of the definition (2.12), we use the spherical coordinates transformation (2.1)-(2.2), and find from the orthogonality (2.8) and (2.11) that which yields (2.14) and ends the proof.
Proposition 2.4. For ą´1 2 and fixed pℓ, q P Υ 8 , we have the following recurrence relations in : and for the GHFs, where " Proof. Recall the three-term recurrence relation of the Laguerre polynomials (cf. [43]): The GHFs with different parameters are connected through the following identity, which finds very useful in the algorithm development.
where the connection coefficients are given by Proof. Recall the property of the generalized Laguerre polynomials (cf. [3] (7.4)): so we can derive the identity from Definition 2.1 and direct calculation.
Remark 2.6. As Γp0q " 8, we can find that in the limiting sense: C " .
Remarkably, for " 0, the GHFs are the eigenfunctions of the harmonic oscillator:´∆`| | 2 , are essential for the error analysis to be conducted in the forthcoming section.

Adjoint generalized Hermite functions in R
Definition 2.8. For ą´1 2 and pℓ, q P Υ 8 , the -dimensional adjoint GHFs are defined by q , ,ℓ p q " where the coefficients t 0 C u are given by (2.21).
Remark 2.9. In light of the connection relation in Proposition 2.5, it is evident that q , ,ℓ p q can be expressed as a linear combination of the counterparts p , ,ℓ p q It is seen from (2.20) (with " 0) that the GHF can be represented as which differs from its adjoint in (2.25) by the sign of the coefficients. Notably, such a subtlety results in an intimate relation between this adjoint pair through the Fourier transform: Moreover, the use of A-GHFs as basis functions in a spectral-Galerkin framework can diagonalise the nonlocal integral fractional Laplacian p´∆q for ą 0. Recall that for ą 0, the fractional Laplacian of P S pR q (functions of Schwartz class) can be naturally defined via the Fourier transform: For 0 ă ă 1, the fractional Laplacian can be equivalently defined by the point-wise formula (cf. [13]): where "p.v." stands for the principle value.

30)
and for ą 0, Moreover, the adjoint GHFs are orthonormal in the sense that for ą 0, Proof. We first show that the GHFs with " 0 are eigenfunctions of the Fourier transform, namely, According to Lemma 9.10.2 of [2], we have that for ą 0, where p q is the Bessel functions of the first kind of order . Then using Definition 2.1 with " 0, and (2.34) with " and " | |, leads to Recall the integral identity of the generalised Laguerre polynomials (cf. [17], p. 820): for ą´1, Thus, taking "`´2 2 in (2.36), we can work out the integral in (2.35) and then obtain from (2.13) with " 0 that This yields (2.32) and ends the proof.
Note that the GHFs are real-valued, so we infer from (2.27) readily that F r p , ,ℓ p qsp q " F´1r p , ,ℓ p qsp q (˚. Thus, we find from (2.30) immediately the following "reversed" form of (2.30).
Corollary 2.11. For ą´1 2 , pℓ, q P Υ 8 and P N 0 , we have F r p , ,ℓ sp q " i`2 q , ,ℓ p q, F´1r q , ,ℓ sp q " p´iq`2 p , ,ℓ p q. (2.39) Remark 2.12. The fractional Sobolev orthogonality (2.32) has profound implications even for the integralorder Laplacian p´∆q with P N. For example, we find from (2.25) with " 1 that the A-GHFs read which are orthogonal with respect to p∇¨, ∇¨q R . However, this attractive property is not valid for the usual Hermite-based methods based on tensorial Hermite functions ś "1 p p q. Thus, it is advantageous to use the A-GHFs for usual Laplacian and bi-harmonic Laplacian (using the A-GHFs with " 2) in R .

Differences and connections with some existing generalisations
There have been some existing generalisations of the usual Hermite polynomials/functions in different senses, so we feel compelled to point out the differences and connections between the GHPs/GHFs herein and some most relevant ones in literature.

(2.44)
It is known that the usual Hermite functions t p u are the eigenfunctions of the Fourier transform. However, this property cannot carry over to the GHFs with " 0. In ( [30] (2.34)), the Fourier transform of p p q p q is expressed in terms of the Kummer hypergeometric function 1 1 p¨q. Here we find from Corollary 2.11 with " 1 the following more informative representation where the one-dimensional adjoint GHFs are given by We remark that the formulation (2.46) requires some simple calculation from (2.25) and (2.40).

2D GHFs versus generalised Hermite bases for Bose-Einstein condensates in [4]
For " 2, the dimensionality of the space ℋ 2 in (2.6) is 2 " 2´0, with the orthogonal basis given by the real and imaginary parts of p 1`i 2 q . In polar coordinates, we have Then by (2.13), the GHFs can be expressed as Note that similar constructions for the 2D GHFs with " 0 have been explored in the computation of the ground states and dynamics of Bose-Einstein condensation (cf. [4]), governed by the Gross-Pitaevskii equation with an angular momentum rotation term: where the constants , ą 0, is the dimensionless angular momentum rotation speed and "´ip BB q "´iB in polar coordinates. The efficient spectral algorithm therein was built upon the constructive basis t e´2 2 p q p 2 qe i u that could diagonalise the Schrödinger operator:´∆`| | 2 . Similar idea was extended to (2.49) in R 3 in cylindrical coordinates by using the tensor product of the 2D basis and the usual Hermite function in the -direction in [4].
In view of Lemma 2.7, the GHFs with " 0 are eigenfunctions of the harmonic oscillator:´∆`| | 2 , so with a proper scaling, the spectral algorithm leads to a diagonal matrix for the scaled harmonic oscillator: ∆`2| | 2 . As we shall show in the late part, our GHFs with " 0 offer a new and efficient tool for the solutions of PDEs involving a more general Schrödinger operator: p´∆q`| | 2 with P p0, 1s and ą´1{2.

3D GHPs versus Burnett polynomials [8]
For " 3, the dimensionality of ℋ 3 in (2.6) is 3 " 2`1. The orthonormal basis in the spherical coordinates " p sin cos , sin sin , cos q takes the form 1 p q " 1 ? 8 p0,0q pcos q; 2 p q " 2`1 ? psin q p , q pcos q cosp q, 2`1 p q " 2`1 ? psin q p , q pcos q sinp q, 1 ď ď , (2.50) where t p , q u are the Gegenbauer polynomials. Then the 3D GHPs/GHFs in Definition 2.1 read more explicit. In fact, for " 0, the GHPs with a scaling turn out to be the Burnett polynomials, which were first proposed by Burnett [8] as follows ,ℓ p q " where is the normalisation constant so that they are orthogonal in the sense As a result, the Burnett polynomials are mutually orthogonal with respect to the Maxwellian ℳp q " 2¨I t is evident that by (2.12) and (2.14) (with " 3 and " 0),

Differences with Hagedorn wavepackets [18] and some other generalisations [46]
The Hagedorn wavepackets first introduced in [18] are deemed as an effective numerical tool in computing quantum dynamics (see [24] for an up-to-date review). They generalise the usual Hermite functions to several dimensions that allow for flexible localisation in position and momentum. According to [23], the Hagedorn wavepackets are constructed from the complex Gaussian function: centred in position P R and momentum P R , where , P Cˆare complex matrices satisfying certain symplecticity conditions, and 0 ă ! 1 is the semiclassical scaling parameter in the time-dependent Schrödinger equation. Then applying the Hagedorns raising operator " p 1 ,¨¨¨, q-fold to 0 p q leads to the Hagedorn wavepackets, which can be represented as a product of the multivariate tensorial Hermite polynomials p q (with an appropriate scaling) and the initial complex Gaussian function as follows where | | " 1`¨¨¨`a nd ! " 1 !¨¨¨!. We refer to [23,24] for more details. Inspired by the work of Hagedorn, the very recent PhD dissertation [46] discussed the extension of the tensorial (usual) Hermite polynomials to the generalised anisotropic Hermite functions of the form where , P Rˆare arbitrary invertible matrices, ą 0 is a parameter and p q " 1 p 1 q¨¨¨p q is a tensor product of the univariate (usual) Hermite polynomials. We refer to [46] for interesting applications in the context of quantum dynamics.
It is evident that our nontensorial GHPs/GHFs are very different from the Hagedorn wavepackets and their variances. We also remark that the tensorial GHFs were briefly discussed in ( [14], p. 278) under a general framework with the weight function ℎ 2 p qe´} } 2 (where ℎ is a reflection-invariant weight function).

GHF approximation of the IFL and the Schrödinger equation
In this section, we implement and analyse the GHF-spectral-Galerkin method for PDEs involving integral fractional Laplacian.

GHF-spectral-Galerkin method for a fractional model problem
As an illustrative example, we consider p´∆q p q`p q " p q in R ; p q Ñ 0 as | | Ñ 8, where P p0, 1q, ą 0, P´pR q, and the fractional Laplacian operator is defined in (2.28)-(2.29). Here, the fractional Sobolev space pR q with real is defined as in [13]. A weak formulation of (3.1) is to find P pR q such that p , q "`p´∆q 2 , p´∆q 2˘R`p , q R " p , q R , @ P pR q. From (2.28), we find readily the continuity and coercivity of the bilinear form p¨,¨q. Then we conclude from the standard Lax-Milgram lemma that the problem (3.2) admits a unique solution satisfying } } pR q ď } }´p R q .
We choose the finite dimensional approximation space spanned by the -dimensional GHFs in Definition 2.1 or equivalently by the A-GHFs in Definition 2.8. However, in view of (2.32), it is advantageous to use the latter as the basis functions, so we define :" span q , ,ℓ p q : 0 ď ď , 1 ď ℓ ď , 0 ď 2 ď´, , ℓ, P N 0 Then, the spectral-Galerkin approximation to (3.2) is to find P such that As with the continuous problem (3.2), it has a unique solution P . In the real implementation, we write and arrange the unknown coefficients in the order and likewise for , but with the components˜, ℓ " p , q , ,ℓ q R . The orthogonality (2.32) implies that the stiffness matrix is an identity matrix. Moreover, in view of the orthogonality of the spherical harmonic basis (cf. (2.8)), the corresponding mass matrix is block diagonal as follows where the entries of each diagonal block can be computed by Thus the linear system of (3.4) can be written as Remark 3.1. With the new basis at our disposal, the above method has remarkable advantages over the existing Hermite approaches (cf. [29,45]). Although the usual one-dimensional Hermite functions are eigenfunctions of the Fourier transform, we observe from (2.28) that the factor | | 2 is non-separable and singular, so the use of tensorial Hermite functions leads to a dense stiffness matrix whose entries are difficult to evaluate due to the involved singularity for ě 2.

Error analysis
Applying the first Strang lemma [42] for the standard Galerkin framework (i.e., (3.2) and (3.4)), we obtain immediately that }´} pR q ď inf P }´} pR q . (3.10) To obtain optimal error estimates, we have to resort to some intermediate approximation results related to certain orthogonal projection. To this end, we consider the 2 -orthogonal projection : 2 pR q Ñ such that p´, q R " 0, @ P . (3.11) From Definition 2.8 and with a change of basis functions, we find readily that :" span p 0, ,ℓ p q : 0 ď ď , 1 ď ℓ ď , 0 ď 2 ď´, , ℓ, P N 0 Thus, we can equivalently write Based on (2.23), we introduce the function space ℬ pR q equipped with the norm where integer ě 0, and`∇ and´∇ are the lowering and raising operators, respectively. The main approximation result is stated below.
Theorem 3.2. Let P p0, 1q. For any P ℬ pR q with integer ě 1, we have Proof. (i). We first estimate the 2 -error. For " 2`1, a direct calculation gives Thanks to the orthogonality (2.15), (2.23)-(3.14) and (3.16), we have that for any ě 0, Then, we derive from (3.13) and (3.17) that If " 2 , we find from (3.14) that (3.16) simply becomes so we can follow the same lines as above to derive the 2 -estimate.
Taking " in (3.10) and using Theorem 3.2, we immediately obtain the following error estimate.

Numerical results
We conclude this section with some numerical results. For the convenience of implementation, we fix the degree of the numerical solution in both radial and angular direction in (3.5), so the numerical solution takes the form , p q " Here, we focus on " 2, 3.
Example 3.5. [Problem (3.1) with a source term.] We next consider (3.1) with the following source functions: The exact solutions are unknown, and we use the numerical solution with " 80, " 20 as the reference solution. For " 2, 3, we plot the maximum errors, in log-log scale, for (3.1) against various in Figure 2(c)-(f), which we take " 0.3, 0.5, 0.7 and fix " 10. As shown in [40], the solution of (3.1) decays algebraically, even for exponentially decaying source terms. Indeed, we observe an algebraic order of convergence.

GHF-spectral-Galerkin method for fractional Schrödinger equations
As a second example, we consider the fractional Schrödinger equation: where P p0, 1s, ą´1{2, the constant ą 0, and the function 0 is given. Here, we focus on the linear equation. Indeed, using a suitable time-splitting scheme, one only needs to solve a linear Schrödinger equation at each time step for some typical nonlinear cases (see, e.g., [4]). We remark that the fractional Schrödinger Equation (3.26) is the model of interest in the study of fractional quantum mechanics, see [22,47], where in [22], this fractional Hamiltonian appeared more reasonable to study the problem of quarkonium. To solve (3.26) efficiently, we adopt the A-GHFs spectral method in space and the Crank-Nicolson scheme in time discretization. Let ∆ be the time-stepping size, and p q « p , ∆ q. Then we look for`1 P pR 2 q such that i´`1∆ ,¯R 2 " where`1 2 " p`1`q{2. We can implement the GHF-spectral scheme as with the problem (3.4), but only need to evaluate the matrix associated with the potential | | 2 . It is a block diagonal matrix To test the accuracy of the proposed method, we add an external source term p , q so that the exact solution is p , q " e´| | 2´. In Figure 3 (a), we plot the maximum errors versus ∆ at " 1, and the second-order convergence is observed. Here we take " 1, " 10, " 50 and different , . We choose the time stepping size to be small so that the error is dominated by the spatial error. In Figure 3(b), we plot maximum errors in the semi-log scale versus various , for which we take " 10, " 1 and different , . We observe that the spatial errors decay exponentially as increases.
Next, we investigate the dynamics of beam propagations as in [47] (where the case " 1 was considered). We take the following incident Gaussian beam as the initial condition: where the constants and are the beam width and the linear chirp coefficient, respectively. In the test, we take " " 1. In Figure 4, we depict the profiles of the real part of the numerical solutions for various , at " 2. Figure 4 (a) shows the solution profile of the usual case with a harmonic potential:´∆`| | 2 for comparison. We observe from the other profiles that the solutions have different peak intensities and singular behaviours, from which we find the smaller the value of , and the stronger the singularity. In fact, some similar observations was made in [47] for the case with " 1.

Müntz-type GHFs with applications to Schrödinger eigenvalue problems
In this section, we introduce the second family of generalised Hermite functions for efficient and spectrally accurate solutions of the Schrödinger eigenvalue problem: where the potential function p q " | | 2 with , being given constants. It is known that (i) if ą´1, all eigenvalues of (4.1) are distinct; (ii) if "´1 or " 0, the spectrum of the Schrödinger operator´1 2 ∆`| | 2 is a continuous one (cf. [15]).
The variational form of (4.1) is to find P R and P 1 pR qzt0u such that ℬp , q :" 1 2 p∇ , ∇ q R`p | | 2 , q R " p , q R , @ P 1 pR q. Numerical solution of (4.12) poses at least two challenges (i) nonpositive definiteness of the variational form and (ii) the singularity of the Coulomb potential. To overcome these, we shall propose an efficient and accurate spectral method by using the Müntz-type GHFs with a suitable parameter " 1 2 , in light of the Coulomb potential.
Define the approximation space , " span p ℋ 1 2 , ,ℓ p q : 0 ď ď , 1 ď ℓ ď 2`1, 0 ď ď , , ℓ, P N 0 ( , where a scaling factor ą 0 is used to enhance the performance of the spectral approximation as in usual Hermite spectral methods in one dimension (see, e.g., [38,44]). The spectral approximation scheme for (4.2) is to find , P R and , P , zt0u such that ℬp , , , q " , p , , , q R 3 , @ , P , . (4.13) In real implementation, we write , p q " "`ˆ0 ,ℓ ,ˆ1 ,ℓ , . . . ,ˆ, ℓ˘, "`ˆ0 1 ,ˆ1 1 ,ˆ1 2 ,ˆ1 3 ,¨¨¨,ˆ1 ,ˆ2 ,¨¨¨,ˆ2`1˘. (4.14) With this ordering, we denote the stiffness and the mass matrices by and , respectively, with the entries given by  Owing to (4.5) and (4.11) with " 1 2 , both the stiffness matrix and the mass matrix are tridiagonal. Consequently, the scheme (4.13) has an equivalent form in the following algebraic eigen-system: " . Interestingly, the matrix`2 8 is diagonal, so we can rewrite (4.15) aś`2 8¯"´`2 8¯, which leads a more efficient implementation. In Figure 5, we plot the errors between the first 30 (counted by multiplicity) smallest numerical eigenvalues and exact eigenvalues in (4.7) versus for fixed " 16 and two different scaling factors (so that the error of the truncation in angular directions is negligible). Observe that the errors decay exponentially in terms of the cut-off number in the radial direction, along which the eigenfunctions are singular. We also see that the scaling parameter affects the convergence rate as the usual Hermite method (cf. [44]).   The corresponding algebraic eigen-system of (4.18) is " . (4.20) In view of orthogonality (4.8) and (4.9), we find that for any P N 0 ,