Exploring Multiscale Quantum Media: High-Precision Efficient Numerical Solution of the Fractional Schrödinger equation, Eigenfunctions with Physical Potentials, and Fractionally-Enhanced Quantum Tunneling

In recent decades, the field of fractional differential equations (FDEs) has developed rapidly due to their wide applicability in physics, engineering, economics, and biology [1, 2, 3, 4, 5, 6, 7, 8, 9]. The best known example of these is the anomalous diffusion equation, which has been experimentally verified to accurately describe diffusion on short time scales [10, 11]. An outstanding question with FDEs is how to properly simulate them with well-converged and understood numerical methods so as to explore and predict such experimental scenarios. For example, anomalous diffusion FDE simulators are used to predict diffusion of pollution in groundwater systems, where the fractional derivatives create a greatly increased probability of pollutants propagating to large distances [12]. Mathematically, this appears as a power-law probability distribution in place of the usual normal or Gaussian one, resulting in the distribution growing faster than the usual $t^{1/2}$ [13]. The anomalous diffusion FDE is a particular example of general reaction-diffusion FDEs which largely remain unsolved and unexplored. In this Article, we focus on the fractional Schrödinger equation (FSE), connected to anomalous diffusion by a complex rotation of the time derivative from the real to the imaginary axis, just as for the integer-order case.

Turning now to possible mathematical solution methods for FDEs, we first note that linear FDE methods such as the Laplace transform, Fourier transform, Mellin transform, Adomain decomposition, and separation of variables may be used effectively in only niche examples [14, 15]. Unfortunately, in general the exact solutions to FDEs are intractable, due in part to the nonlocal nature of fractional derivatives. Such nonlocality can be observed in the many definitions of fractional derivatives. For example the Riesz fractional derivative is a particularly useful form for physics readers, since it generalizes from the Fourier transform $\partial/\partial x\to ik$ as follows:

$$(\partial^{\alpha}/\partial x^{\alpha})f(x)\equiv\partial_{x}^{\alpha}f(x)\equiv\int dk\,e^{ikx}(-|k|^{\alpha})\tilde{f}(k)\,.$$

(1)

In effect, this nonlocality creates an integrodifferential equation. Thus, numerical methods are a vital component to understanding FDEs. Direct implementation of fractional derivatives to discrete numerical schemes is subject to creating errors and inaccuracies just as with integer derivatives. For example, famously the obvious direct discretization of the usual integer-order Schrödinger equation via finite-differencing does not conserve the normalization [16]. Many FDE methods developed so far, while useful for back of the envelope brainstorming, suffer from such issues. Thus, in order to make accurate and reliable predictions for multiscale experiments, rigorous and careful development of such numerical methods is desirable. The methods we develop in this Article will be generally useful toward this goal. As a clear demonstration and for context we will focus on solving the eigenfunction problem for the FSE in a variety of physically relevant potentials. Eigenfunctions provide a complete basis for all solutions for linear FDEs just like integer order linear differential equations.

One tracing of the historical origin of the FSE can be pieced together as follows. The original Schrödinger equation was reinterpreted as a path integral over Brownian motion by Feynman [17]. In Brownian motion, the paths can be written as non-differentiable, self-similar curves with fractional dimensions that may differ from the Cartesian embedding in e.g. 3D. Laskin then constructed the FSE as a path integral over Lévy flights [18]. This construction is the natural generalization over paths with fractional dimension differing from Brownian motion. This change in the dispersive properties of the underlying equation creates derivatives of fractional order. The particular type of derivative that arises from this construction is the Riesz fractional derivative, Eq. (1). Although there are many useful definitions of fractional derivatives, such as the Caputo derivative used in some anomalous diffusion FDE simulation methods [19], we will focus on the Riesz derivative for the purposes of this Article.

In the context of Laskin’s many seminal contributions to the formulation and understanding of the FSE [18, 20, 21, 22] we point out that he naturally used a common approach to solving the integer-order Schrödinger equation $\alpha=2$ for piece-wise potentials, namely matching the wavefunction and its first derivative at step-function boundaries. This is a general method for linear ODEs commonly taught for a wide range of reaction-diffusion type equations in which the number of required matching conditions at each step-function is equal to one less than the order of the highest derivative. Although not a universal method it applies in a wide range of cases. The nonlocal nature of fractional derivatives is one definite exception to this approach, since regions bounding the step functions cannot be treated independently [23]. This led to a misunderstanding of evanescent waves in Laskin’s work [18] which we shall correct. The essential elements of our numerical approach can be applied to a variety of reaction-diffusion and other ODE and even PDE contexts addressing decay of a solution in a barrier. As particular applications to the FSE we will consider first the particle on a ring problem showing that the functional form is in this case independent of the order of the derivative, while the eigenvalue does depend on $\alpha$. We then consider the quantum harmonic oscillator, where we show that decay into the barrier is greatly slowed for non-integer orders, enhancing the evanescent waves. Finally for the finite well and double well we show this enhanced evanescent wave effect naturally leads to higher tunneling rates due to larger penetration into the barrier, which we call fractionally-enhanced quantum tunneling, an effect one expects to observe in multiscale media in quantum systems.

The one-dimensional FSE takes the following form:

$$i\hbar\partial_{t}\Psi=-C_{\alpha}\frac{\hbar^{2}}{2m}\partial_{x}^{\alpha}\Psi+V(x)\Psi\,,$$

(2)

where the Reisz fractional derivative of order $\alpha$ may be here defined more succinctly than Eq. (1) in terms of the Fourier transform $F$

$$\partial_{x}^{\alpha}f(x)=\mathcal{F}^{-1}(-|k|^{\alpha}\mathcal{F}(f(x)))\,,$$

(3)

where $C_{\alpha}$ is a constant with physical units that correct those imposed by the fractional derivative. Clearly, the Reisz fractional derivative is a Fourier multiplier type operator and given a fractional order of $\alpha=2$ we recover the ordinary second order derivative.

To simplify our analysis we rescale our units such that the FSE simplifies to the following.

$$i\partial_{t}\Psi=-\frac{1}{2}\partial_{x}^{\alpha}\Psi+V(x)\Psi=\hat{H}\Psi\,.$$

(4)

Given this construction we may define a time evolution operator by integrating out the first derivative with respect to time,

$$\Psi(x,\Delta t)=e^{-i\hat{H}\Delta t}\Psi(x,0)\,.$$

(5)

The eigenfunction problem may then be defined as those functions for which application of our Hamiltonian applies only a scaling factor to our function:

$$\hat{H}\phi_{n}(x)=E_{n}\phi_{n}(x)\,.$$

(6)

The set of linearly independent eigenfunctions which fulfills this property defines a basis from which any solution may be decomposed [24]. The goal of this paper is to construct an effective means by which we may numerically construct eigenfunctions within this set. The means by which we will construct these eigenfunctions is via imaginary time evolution through the Fourier split step method, as discussed in Sec. 4.

Before discussing the Fourier split step method, a special function with many connections to fractional derivatives must be briefly discussed. The Mittag-Leffler function arises naturally in the solution of many fractional order integral equations or fractional order differential equations and has diverse application in the study of fractional relaxation and fractional diffusion [25, 26]. Although the eigenfunctions generated to solve the fractional Schrodinger equation are not quantitatively described by the Mittag-Leffler function, certain behavior attributed to the Mittag-Leffler function and its family of functions appears when examining the eigenstates, especially the decay patterns in the classically forbidden regions of the eigenfunctions.

The Mittag-Leffler function can be understood as a type of generalization to the exponential function and is defined with one or two parameters via the following sums:

$$E_{q}(x)=\sum_{k=0}^{\infty}\frac{x^{k}}{\Gamma(qk+1)}\,,$$

(7)

$$E_{q,\beta}(x)=\sum_{k=0}^{\infty}\frac{x^{k}}{\Gamma(qk+\beta)}\,.$$

(8)

Where the exponential function is the eigenfunction of the derivative operator, the Mittag-Leffler function is the eigenfunction of the Caputo fractional derivative, a fractional derivative connected to the fractional Rietz operator we study in this article [27]. As a practical example of the type of behavior that the Mittag-Leffler function portrays, Fig. 1 compares a Gaussian and its Mittag-Leffler analog.

Clearly the Gaussian function is simply the Mittag-Leffler function given that $q=1$ with the parameter $-x^{2}$. Additionally, close to the origin the functions are similar with slight divergences of slower or faster decay. However, the important piece is in the tails of Fig. 1(b). Away from the central curve, the tails present a stark difference for non-integer $q$. Examining the case of $q=0.9$, they decay significantly slower than the Gaussian case and in fact the decay appears to slow down as the tail extends further. Furthermore, for $q=1.1$ we observe the main difference is the appearance of a node before it again slowly decays towards zero. In both cases we find the tails decay significantly slower for non-integer $\alpha$ as compared to the integer case.

To summarize, where typically we would observe exponentially decaying tails in integer order eigenfunctions, what we find in our generated eigenfunctions with the FSE is Mittag-Leffler-like decaying tails. Although not exactly the same, certain stark features very reminiscent of Mittag-Leffler functions and related functions may be observed.

Recall that imaginary time evolution takes the form

$$f_{0}(x,-it)=e^{-\hat{H}t}f_{0}(x)=\sum_{n}a_{n}e^{-E_{n}t}\psi_{n}(x)$$

(9)

with $\psi_{n}(x)$ the eigenfunctions of a Hamiltonian and $a_{n}$ the initial weights. All higher energy components decay exponentially more rapidly in imaginary time, resulting in the ground state. To obtain higher eigenfunctions, one simply subtracts lower eigenfunctions at every renormalization step. A modification to this method must be used to treat nearly degenerate or degenerate eigenstates, as one encounters e.g. in the double well potential. For the double well it is sufficient to constrain the parity to be even in order to obtain the ground state, and use alternating parity constraints for excited states.

The difficulty then is not how to generate the eigenfunctions, but how to evolve in time given a fractional derivative in the Hamiltonian. Although the Reisz fractional derivative is quite complicated and difficult to apply in position space, its application in momentum space is straightforward as per Eq. (1). Thus we have a portion of our Hamiltonian which is easy to apply in position space and difficult to apply in momentum space, the potential; and a portion which is easy to apply in momentum space and difficult to apply in position space, the fractional derivative. The Fourier split step method is a well-known solution to this circumstance, which we here write out plainly for the fractional case.

Given that we may split our Hamiltonian into these components where $\hat{M}$ is the operator that is easy to apply in the momentum domain, and $\hat{P}$ is the operator that is easy to apply in the position domain, the problem is that we may not immediately separate these operators within the exponential function as these operators in general do not commute. However, by ignoring this and splitting them regardless we may construct an approximation that is 1st order in terms of the time step, $\Delta t$,

$$\hat{H}=\hat{M}+\hat{P}\,,$$

(10)

$$e^{-\hat{H}\Delta t}=e^{-\hat{M}\Delta t}e^{-\hat{P}\Delta t}+\mathcal{O}(\Delta t)\,.$$

(11)

As per the usual split step method, we Fourier transform the state before applying $\hat{M}$, then inverse Fourier transform, apply $\hat{P}$, and obtain the first order approximation. Given that our operators have the following form the complete procedure is as follows:

$\displaystyle\hat{P}=V(x)$

$\displaystyle\hat{M}=\frac{1}{2}|k|^{\alpha}\,,$

(12)

$$\psi(x,-i\Delta t)=\mathcal{F}^{-1}\left(e^{-\frac{1}{2}|k|^{\alpha}\Delta t}\mathcal{F}\left(e^{-V(x)\Delta t}\psi(x,0)\right)\right)+\mathcal{O}(\Delta t)\,.$$

(13)

We may then increase the accuracy of this scheme by adding more steps and changing the constants on the operators such that the exponential is applied more accurately. A general construction of a splitting scheme is

$$e^{\hat{A}+\hat{B}}\approx\prod_{k=1}^{n}e^{a_{k}\hat{A}}e^{b_{k}\hat{B}}\,.$$

(14)

An example of a different scheme that may be recovered via this procedure is the second order Suzuki-Trotter approximation, which is also a second order approximation in time:

$$e^{(\hat{A}+\hat{B})\Delta t}=e^{\frac{1}{2}\hat{A}\Delta t}e^{\hat{B}\Delta t}e^{\frac{1}{2}\hat{A}\Delta t}+\mathcal{O}(\Delta t^{2})\,.$$

(15)

However, there is a significant subtlety here – the condition on which constants $a_{n}$ and $b_{n}$ are allowed in order to keep the scheme unconditionally stable for the case of imaginary time evolution. All constants $a_{n}$ and $b_{n}$ must be positive for all $n$. Unlike the real time evolution case where the applied complex exponential functions are rotations in complex space, imaginary time evolution has either a decaying or growing exponential. If $a_{n}$ or $b_{n}$ is negative then higher energy eigenfunctions grow significantly faster during specific steps than the lower energy eigenfunctions that we wish to isolate. Furthermore, machine error will produce contributions of very high order eigenfunctions that will grow significantly faster and lead to numerical instability. This problem is the same as that of evolving the diffusion equation backwards in time. The problem that then arises is that splitting methods with approximation orders greater than 2nd order must necessarily have negative splitting coefficients [28]. We can fix this is by allowing our splitting coefficients to not only be real valued but also complex valued. Given this adjustment it is easy to show that higher-order splitting schemes become available at the cost of the additional overhead of using complex numbers. Given this additional freedom, we choose to use a sixth order scheme that takes 8 steps. The coefficients from Eq. (14) are found in Table 1.

To summarize, our FSE code is based on a sixth-order split-step method within the framework of the Suzuki-Trotter approximation. This approach facilitates the efficient application of the Riesz fractional derivative by leveraging the Fast Fourier Transform (FFT). A notable aspect of using FFT in our simulations is the implicit assumption of periodic boundary conditions. However, for scenarios involving bound states, open boundary conditions can be effectively emulated by simply expanding the system size. This expansion ensures that the boundary effects do not significantly influence the dynamics within the region of interest.

In this context, when evolving the Schrödinger equation in imaginary time, which corresponds to the real-time evolution of a diffusion equation with an added potential term, our approach employs the Fourier split-step method. This method, by its very nature, is unconditionally stable when applied to linear, time-independent problems. Its stability stems from the precise handling of the linear part of the equation in Fourier space, complemented by the treatment of the nonlinear or potential-dependent part in real space within each discrete time step. This characteristic of unconditional stability is particularly advantageous for simulations aiming to reach the ground state, as does not induce a restriction of the values of $\Delta t$ and $\Delta x$, this then permits the use of larger time steps without risking numerical instability, thereby facilitating efficient convergence over extended imaginary time scales.

However, the integration of certain types of potentials into the simulation imposes additional considerations related to spatial and frequency resolutions. Specifically, potentials characterized by large magnitudes, abrupt spatial variations, or significant nonlinearities necessitate a meticulous choice of spatial grid size ($\Delta x$) and time step ($\Delta t$). The presence of steep potential gradients calls for a denser spatial grid to accurately capture these variations, while high-frequency components induced by the potential require careful management in the frequency domain to prevent aliasing effects. Therefore, while the primary stability concern of the Fourier split-step method in the temporal domain is addressed, the introduction of complicated potentials demands a nuanced approach in spatial discretization and frequency domain analysis. This is essential to ensure the accuracy and reliability of the simulation results, particularly in quantum systems where the potential landscape critically influences the system dynamics.

In this study, an extensive convergence analysis was conducted to identify the optimal grid sizes that effectively minimize finite grid errors below the threshold relevant to our function analyses. For the entirety of the results presented in this paper, we standardized the grid size and time step at $\Delta x=10^{-2}$ and $\Delta t=10^{-2}$. An exception is made in the case of the finite well problem, where an additional step is incorporated to accurately address the sharp potential variations. Details of such adjustments are elaborated in the corresponding sections of the paper.

In line with our commitment to open science and reproducibility, the complete source code for our Fourier Split-Step simulations, as described in this paper, is made publicly available. Interested researchers can access, review, and utilize the code through our GitHub repository, titled ”FractionalEigenstate.” This repository includes detailed documentation and example scripts, providing a comprehensive guide for users to adapt and apply the code to their specific research needs. For access to the repository, please refer to our citation below [29].

To examine the effectiveness of this algorithm on fractional order quantum systems we first investigate one of the very few systems in which the analytical energies and eigenstates are well known, the particle on a ring system. This will enable a careful evaluation of our code’s precision.

We define the particle on a ring to be a constant potential, here chosen to be zero, under periodic boundary conditions, with domain $x\in[-1,1]$. The Hamiltonian is thus

$$\hat{H}=-\frac{1}{2}\partial_{x}^{\alpha}.$$

(16)

Given that the Reisz fractional derivative is a Fourier multiplier similar to integer derivatives, the eigenfunctions of such an operator by itself are plane waves. As we cast our eigenfunctions to the real plane to ensure uniqueness these eigenfunctions are then just sines and cosines similar to the particle on a ring problem with the second derivative. Because of this, the eigenfunctions of this system are invariant of fractional order. The only aspect of the eigen-problem that does change is the energy, which we can obtain by imaginary time evolution, as per Sec. 4, as well as analytically for this special case:

$$\partial_{x}^{\alpha}\sin(kx)=-|k|^{\alpha}\sin(kx)\,,$$

(17)

$$E_{n}=\frac{1}{2}\left(\pi\lceil n/2\rceil\right)^{\alpha}\,,$$

(18)

where $\lceil*\rceil$ denotes the ceiling function. This is due to the degeneracy of the energies within the system. The reader will note that the usual integer order expression for the energy is recovered when $\alpha=2$. The eigenstates of the system are given by

$$\psi_{n}(x)=\begin{cases}\cos{\pi nx}&\text{n even}\\ \sin{\pi nx}&\text{n odd}\,,\end{cases}$$

(19)

independent of fractional order. To ascertain the effectiveness of this algorithm, we first compare the generated first excited state for the particle on a ring system to its analytical eigenstate across a few fractional orders. For this examination, we define our grid in the domain [-1,1] with 480 evenly spaced points, a choice strategically made to enable the representation of the 39th excited state within the constraints of the Nyquist sampling theorem. The results of this analysis are presented in Figure 2. Interestingly higher order systems tend to converge faster and yield more accurate results than their lower order counterparts. Although not universally applicable, in many systems, an increase in the fractional order leads to a greater energy spacing between eigenstates. This increase in energy spacing inherently accelerates the convergence process via imaginary time evolution.

We now focus on evaluating the maximum error in the pointwise approximation of each eigenstate, with a particular interest in how this error varies in relation to the ordinal number of the eigenstate. To conduct this analysis, we opt for a discretization of the domain into 480 uniformly distributed points. The rationale behind choosing this specific number of points is grounded in the principles of the Nyquist sampling theorem. According to this theorem, to accurately represent a signal with a frequency of $f=\pi n$, a minimum of $4\pi n+1$ points is required. Hence, for effectively capturing the nuances of the 39th excited state, which is represented by the function $\cos{39\pi x}$, a minimum of 479 points is necessary. Therefore, our selection of 480 points ensures that we adequately sample the highest frequency state in our analysis, minimizing the potential for aliasing errors.

Finally we may compare the calculated energies from our generated eigenstates to the analytical energies for this system.

Having established the efficacy of our code for the eigenvalue problem, we proceed to treat three classic quantum problems of physical interest for fractional systems: the quantum harmonic oscillator, the finite well, and the double well. We show how the new parameter $\alpha$ in the FSE creates new features. The double well in particular contains the essential elements of quantum tunneling between discrete energy levels, and has some quite surprising new aspects, which we call fractionally-enhanced quantum tunneling.

In this article we have constructed an algorithm capable of solving for the eigensystem of the fractional Schrödinger equation with an arbitrary physical potential up to machine precision. Given this tool we have then examined a few physically important systems and discovered that the exponential evanescent waves distinctive to classically forbidden regions are replaced with many of the same properties of the fractional generalization to the exponential function, the Mittag-Leffler function. We showed the fractional order of the system can be used to control and enhance tunnelling via the non-locality and enhanced evanescent waves associated with the fractional Schrödinger equation. We have termed this effect fractionally-enhanced quantum tunneling. This algorithm is a tool that may then be used by experimentalists and applied mathematicians to further investigate the properties of multi-scale materials and discover unique and novel ways to design quantum devices and materials.

Acknowledgements

We acknowledge useful discussions with Joel Been and Mingzhong Wu. This work was supported in part by the NSF under grant number DMR-2002980.