# Fractional Schrödinger Equation and Time Dependent Potentials

E. C. Gabrick${}^{1},$111Corresponding Author
e-mail<EMAIL_ADDRESS>E. Sayari1, A. S. M. de Castro1,2, J. Trobia3, A.
M. Batista1,3, E. K. Lenzi1,2 1Graduate Program in Science, State University
of Ponta Grossa, Ponta Grossa, PR, Brazil
2Department of Physics, State University of Ponta Grossa, Ponta Grossa, PR,
Brazil
3Department of Mathematics and Statistics, State University of Ponta Grossa,
Ponta Grossa, PR, Brazil

###### Abstract

We investigate the solutions for a time-dependent potential by considering two
scenarios for the fractional Schrödinger equation. The first scenario analyzes
the influence of the time-dependent potential in the absence of the kinetic
term. We obtain analytical and numerical solutions for this case by
considering the Caputo fractional time derivative, which extends Rabi’s model.
In the second scenario, we incorporate the kinetic term in the Schrödinger
equation and consider fractional spatial derivatives. For this case, we
analyze the spreading of the Gaussian wave package under the action of the
time and spatial fractional differential operators.

anomalous spreading, fractional dynamics, fractional quantum mechanics

##  I Introduction

The meaning of the operator $d^{\nu}y/dx^{\nu}$ with $\nu$ integer is well
known and has a profound physical background [2]. The challenge is to
understand what this operator means when $\nu$ is any number (positive or
negative, real or complex) [3] or even a function [4]. This problem can be
dated from a letter of L’Hôpital to Leibniz in 1695, where he asked him what
the operator $d^{\nu}y/dx^{\nu}$ is when $\nu=1/2$ [3]. Since then, many
researchers have dedicated themselves to this problem, for example, Euler,
Lagrange, Laplace, Fourier, and others [2], giving rise to the fractional
calculus [5].

Nowadays, fractional calculus has quickly become a new efficient mathematical
tool to analyze different properties of systems, in general, by extending the
differential operators by incorporating non-integer indexes and, in
particular, connecting them with experimental results [6, 7, 8, 9]. In this
manner, it is possible to investigate many situations with a simple extension
that may incorporate memory effects, long-range correlations, and other
effects in complex systems [10]. For instance, in complex viscoelastic media
[11, 12], electrical spectroscopy impedance [13, 14, 15], wave propagation in
porous media [16, 17], microflows of viscoelastic fluids [18], and gas
transport in heterogeneous media [19, 20, 21]. It has also been used in other
physics branches to extend several partial differential equations to cover and
bring new possibilities for applications in different scenarios [2]. One of
them is quantum mechanics, which has been extended by incorporating spatial
and time fractional differential operators [22, 23]. In this context, we have
the pioneers works of N. Laskin [24, 25, 26], which lead us to a fractional
Schrödinger equation, and that has been followed by other extensions
incorporating fractional differential operators in time, and space [3] as well
as non-local terms [27, 28] and constraints among the different spatial
coordinates (comb-model) [29, 30, 31]. These extensions of the Schrödinger
equation have also been analyzed by considering different choices of
potential, such as delta potentials [32], constant or linear potentials [33]
and for some time dependent potentials [34]. It is worth mentioning that, from
the analytical and numerical point of view, it is a challenge to obtain
solutions when fractional time derivatives are considered.

Our goal in this work is to investigate the implications of considering time
dependent potentials in the following fractional Schrödinger equation [35]

$\displaystyle i^{\alpha}\hbar_{\alpha}\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi(\vec{r},t)$ $\displaystyle=$
$\displaystyle\widehat{H}(t)\psi(\vec{r},t)\;,$ (1)

where the fractional differential operator is the Caputo fractional time
derivative, defined as follows [3]:

$\displaystyle\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi(\vec{r},t)=\frac{1}{\Gamma\left(1-\alpha\right)}\int_{0}^{t}dt^{\prime}\frac{1}{(t-t^{\prime})^{\alpha}}\frac{\partial}{\partial
t}\psi(\vec{r},t),$ (2)

for $0<\alpha<1$. We employ analytical and numerical approaches to analyze Eq.
(1). For the last one, we consider the ﬁnite difference method [66, 67, 68].
It should be mentioned, as discussed in Ref. [35], that we can also extend the
Schrödinger equation as follows:

$\displaystyle i\hbar_{\alpha}\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi(\vec{r},t)$ $\displaystyle=$
$\displaystyle\widehat{H}(t)\psi(\vec{r},t).$ (3)

Equations (1) and (3) are two possible extensions of the Schrödinger equation.
However, when performing a Wick rotation, the imaginary unit is raised to the
same power as the time coordinates for Eq. (1). Another point between the two
equations involves the temporal behavior of the solution, which for the first
case, is more suitable than the second one that decreases or grows with time
instead of a sinusoidal behavior. For these reasons point out in Ref. [35], we
consider Eq. (1) in our developments. It is also interesting to mention the
similar appearance between the Schrödinger and diffusion equations. This
similarity between these equations is a consequence of the stochastic
processes behind these equations, which can be evidenced by Feynman’s path
integral formulation [36], and transformed into Wiener’s path integral, which
is the integral over the path of Brownian motions. It has also motivated
different extensions motivated by other aspects, which include Lévy
distributions [26], comb-model [37, 38], among others. In addition, these
extensions of the Schrödinger equations have been considered in problems
related to optica [39], solutions for free-particle [40], optical solitons
[41] and others [42, 43, 44, 45, 46].

By using Eq. (1), we consider a two-level system with a time dependence on the
potential and restricted to a one-dimensional wave function $\psi(x,t)$,
without any loss of generality, and $\hbar_{\alpha}$ is an arbitrary time
constant used to replace the Planck constant (see Ref. [35] for more details).
As mentioned before, the difference between the definitions given by Eq. (1)
and Eq. (3) is in the imaginary unit. Both equations violate the probability
conservation law [47]. However, the probability related to Eq. (1) may
increase and reach a constant value $1/\alpha^{2}$ as discussed in Ref. [48]
and the probability associated with Eq. (3) decays to zero [47]. It is worth
mentioning that the two-level systems are very interesting because the
simplicity and richness of results [49] have been used to study spin 1/2-like
[50], magnetic resonance [51], quantum computation [52], unitary evolution
[53], and others [54]. In some cases, the two-level systems are analytical
soluble, mostly when the Hamiltonian is unperturbed. However, perturbed
Hamiltonians are particularly interesting, mainly in the presence of an
electromagnetic field [55]. In situations like that, i.e., the time-dependent
Hamiltonian, the exact solution is rare; one famous example is the Rabi
problem [56]. Inspired by the Rabi problem and electromagnetic fields
perturbation, we consider two distinct cases for the Hamiltonian in Eq. (1).
The first one considers

$\widehat{H}(t)=\left(\begin{array}[]{cc}E_{1}&\gamma e^{i\omega t}\\\ \gamma
e^{-i\omega t}&E_{2}\end{array}\right),$ (4)

which corresponds to a two-level system, where $E_{1}$ and $E_{2}$ are the
eigenvalues, and $\gamma$ is the amplitude of the external field with
frequency equals $\omega$. In the second case, we consider the Hamiltonian
given by

$\widehat{H}(t)=\left(\begin{array}[]{cc}\widehat{p}^{\,\mu}/(2m)&\gamma
e^{i\omega t}\\\ \gamma e^{-i\omega
t}&\widehat{p}^{\,\mu}/(2m)\end{array}\right),$ (5)

which incorporates a kinetic term and consequently a spatial dependence in our
problem. Note that the kinetic terms have the power $\mu$, which can be
related to a spatial fractional derivative, i.e.,
${\cal{F}}^{-1}\\{|p|^{\mu}\widetilde{\psi}(p,t)\\}\equiv-\partial^{\mu}_{|x|}\psi(x,t)$,
where ${\cal{F}}\\{\psi(x,t);k\\}=\widetilde{\psi}(k,t)$ (and
${\cal{F}}^{-1}\\{\tilde{\psi}(k,t);x\\}=\psi(x,t)$) corresponds to the
Fourier transform, respectively. This definition corresponds to the Riesz
derivative [57, 58].

Aiming to understand the influence of fractional order in Schrödinger
equation, we made the developments for standard quantum mechanics in Sec. II
and for the fractional operators in Sec. III. We obtain analytical and
numerical solutions for these Hamiltonians and analyze the spreading behavior
of the wave package in different conditions. Finally, we present our
discussions and conclusions in Sec. IV.

##  II Schrödinger Equation

The standard Schrödinger equation is an specific case of Eq. (1) or Eq. (3)
with $\alpha=1$. To understand the effects of $\alpha\neq 1$ in quantum
dynamics, we first analyze the results obtained from the standard case. In
this sense, let us start our analysis by reviewing the results obtained for
the standard Schrödinger equation, i.e.,

$i\hbar\frac{\partial}{\partial
t}\psi(\vec{r},t)=\widehat{H}\psi(\vec{r},t)\;,$ (6)

where $\widehat{H}$ is the Hamiltonian operator, $\psi(\vec{r},t)$ is the wave
function, $i$ is the imaginary unit, and $\hbar$ is the Planck constant [49],
which, for simplicity, we consider $\hbar=1$. Equation (6) will be analyzed
first by considering the Hamiltonian given by Eq. (4) which corresponds to a
two-level system, as previously discussed. Equation (4) has been applied in
several situations, such as a two-level system interacting with light field
[48]. After, we incorporate kinetic terms in Eq. (4) by performing the
following change $E_{1}\rightarrow\widehat{p}^{\,2}/\left(2m\right)$ and
$E_{2}\rightarrow\widehat{p}^{\,2}/\left(2m\right)$, which implies

$\widehat{H}=\left(\begin{array}[]{cc}\widehat{p}^{\,2}/\left(2m\right)&\gamma
e^{i\omega t}\\\ \gamma e^{-i\omega
t}&\widehat{p}^{\,2}/\left(2m\right)\end{array}\right)\;.$ (7)

Equation (7) is equivalent to considering the particular case $\mu=2$ in Eq.
(5), i.e., it considers the kinetics terms with an integer index. After
analyzing the standard Schrödinger equation which emerges from these cases, we
consider the fractional extensions of these cases and analyze the implications
for spreading the wave package, particularly the case $\mu\neq 2$. Equations
(4) and (7) allows us to consider that the wave function has the following
form

$\displaystyle\psi=\left(\begin{array}[]{cc}\psi_{1}\\\
\psi_{2}\end{array}\right)\;,$ (10)

with $\psi_{1}$ and $\psi_{2}$ are obtained by solving the Schrödinger
equation for each case.

Now let us consider the first case correspondent to the Hamiltonian, defined
in terms of Eq. (4) and solutions $\psi_{k}=\psi_{k}(t)$. For the initial
condition, we analyze the situation in which only one state is populated
initially, i.e., the initial condition is given by $\psi_{1}(0)=1$ and
$\psi_{2}(0)=0$. The problem concerns obtaining the probability transition
between states after applying the external field. We find these probabilities
by solving Eq. (6), i.e.,

$\displaystyle i\frac{\partial}{\partial t}\psi_{1}(t)$ $\displaystyle=$
$\displaystyle E_{1}\psi_{1}(t)+\gamma e^{i\omega t}\psi_{2}(t),$ (11)

and

$\displaystyle i\frac{\partial}{\partial t}\psi_{2}(t)$ $\displaystyle=$
$\displaystyle E_{2}\psi_{2}(t)+\gamma e^{-i\omega t}\psi_{1}(t),$ (12)

in which $|\psi_{1}(t)|^{2}+|\psi_{2}(t)|^{2}=1$ is always verified. This case
admits an analytical solution; for example, see Ref. [49] and, in particular,
when $\omega=0$, it is given by

$\displaystyle\psi_{1}(t)=\frac{1}{2}\left(1+\frac{E_{1}-E_{2}}{\sqrt{(E_{1}-E_{2})^{2}+4\gamma^{2}}}\right)e^{-i\gamma_{+}t}+\frac{1}{2}\left(1-\frac{E_{1}-E_{2}}{\sqrt{(E_{1}-E_{2})^{2}+4\gamma^{2}}}\right)e^{-i\gamma_{-}t},$
(13)

and

$\displaystyle\psi_{2}(t)=\frac{\gamma}{\sqrt{(E_{1}-E_{2})^{2}+4\gamma^{2}}}\left(e^{-i\gamma_{+}t}-e^{-i\gamma_{-}t}\right),$
(14)

where
$\gamma_{\pm}=\left(E_{1}+E_{2}\pm\sqrt{(E_{1}-E_{2})^{2}+4\gamma^{2}}\right)/2$.

Figure 1(a) illustrates the result obtained by considering the external field
constant, i.e., $\omega=0$. It is possible to verify that the system has
oscillation between two levels. It is interesting to note in Fig. 1(a) that
for the parameters used, $\psi_{1}$ is predominant over $\psi_{2}$. On the
other hand, when we consider an oscillatory external field, $\omega\neq 0$,
the system oscillates between the two states as shown in Fig. 1(b). The result
for $\omega\neq 0$ is numerically obtained by solving Eqs. (11) and (12). When
$0<\omega<1$, the amplitude of $|\psi_{2}|^{2}$ tends to 1, founding this
value in $\omega=1$. By the other hand, for $\omega>1$ the $|\psi_{2}|^{2}$
value oscillate asymptotically to zero while $|\psi_{1}|^{2}$ oscillate in 1
direction.

Figure 1: Probability of finding the system in $\psi_{1}(t)$ state given by
the red line and in the $\psi_{2}(t)$ by the blue line. The panel (a) is for
$\omega=0$ and (b) for $\omega=1$. We consider $\gamma=1$, $E_{1}=1$ and
$E_{2}=2$.

Similar analysis can be performed by considering the second case, i.e., for
the Hamiltonian defined in terms of Eq. (7). The equations for the wave
functions $\psi_{1(2)}(x,t)$ read as

$\displaystyle i\frac{\partial}{\partial t}\psi_{1}(x,t)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\psi_{1}(x,t)+\gamma e^{i\omega t}\psi_{2}(x,t),$ (15)

and

$\displaystyle i\frac{\partial}{\partial t}\psi_{2}(x,t)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\psi_{2}(x,t)+\gamma e^{-i\omega t}\psi_{1}(x,t),$ (16)

where, for simplicity, we assume $m=1$. The solution for this case can be
found by applying the Fourier transform
($\widetilde{\psi}_{1,2}(k,t)={\cal{F}}\\{\psi_{1,2}(x,t);k\\}$ and
$\psi_{1,2}(x,t)={\cal{F}}^{-1}\\{\widetilde{\psi}_{1,2}(k,t);x\\}$), as
defined before) in Eqs. (15) and (16) yielding

$\displaystyle i\frac{\partial}{\partial t}\widetilde{\psi}_{1}(k,t)$
$\displaystyle=$
$\displaystyle\frac{1}{2}k^{2}\widetilde{\psi}_{1}(k,t)+\gamma e^{i\omega
t}\widetilde{\psi}_{2}(k,t),$ (17)

and

$\displaystyle i\frac{\partial}{\partial t}\widetilde{\psi}_{2}(k,t)$
$\displaystyle=$
$\displaystyle\frac{1}{2}k^{2}\widetilde{\psi}_{2}(k,t)+\gamma e^{-i\omega
t}\widetilde{\psi}_{1}(k,t)\;.$ (18)

By performing some calculations, it is possible to show that the solution
$\widetilde{\psi}_{2}(k,t)$ is related to the solution
$\widetilde{\psi}_{1}(k,t)$ as follows:

$\displaystyle\widetilde{\psi}_{2}(k,t)=-i\gamma\int_{0}^{t}dt^{\prime}e^{-\frac{1}{2}ik^{2}(t-t^{\prime})}e^{-i\omega
t^{\prime}}\widetilde{\psi}_{1}(k,t^{\prime})\;,$ (19)

for which we assume the initial condition $\widetilde{\psi}_{2}(k,0)=0$.
Furthermore, this relation implies that

$\displaystyle i\frac{\partial}{\partial
t}\widetilde{\psi}_{1}(k,t)=\frac{1}{2}k^{2}\widetilde{\psi}_{1}(k,t)-i\gamma^{2}\int_{0}^{t}dt^{\prime}e^{-\frac{1}{2}ik^{2}(t-t^{\prime})}e^{-i\omega(t-t^{\prime})}\widetilde{\psi}_{1}(k,t^{\prime})\;.$
(20)

Note that the last term present in Eq. (20) is a nonlocal term and the kernel
has a nonsingular dependence on the variable $t$. It is worth mentioning that
the nonsingular kernels have been successfully applied in many situations,
such as the ones presented in Refs. [59, 60, 61, 62, 63, 64].

Equation (20) can be solved by using the Laplace transform, yielding

$\displaystyle\widetilde{\psi}_{1}(k,t)=e^{-\frac{1}{2}i\left(k^{2}-\omega\right)t}\left[\cos\left(\frac{1}{2}t\sqrt{\omega^{2}+4\gamma^{2}}\right)-\frac{i\omega}{\sqrt{\omega^{2}+4\gamma^{2}}}\sin\left(\frac{1}{2}t\sqrt{\omega^{2}+4\gamma^{2}}\right)\right]\widetilde{\varphi}_{1}(k)\;,$
(21)

where $\widetilde{\psi}_{1}(k,0)=\widetilde{\varphi}_{1}(k)$ is the initial
condition for $\psi_{1}(x,t)$. Applying the inverse of the Fourier transform,
we obtain that

$\displaystyle\psi_{1}(x,t)=\Xi_{1}(t)\int_{-\infty}^{\infty}dx^{\prime}{\cal{G}}(x-x^{\prime},t)\varphi_{1}(x^{\prime}),$
(22)

where

$\displaystyle\Xi_{1}(t)=e^{\frac{i}{2}\omega
t}\left[\cos\left(\frac{1}{2}t\sqrt{\omega^{2}+4\gamma^{2}}\right)-\frac{i\omega}{\sqrt{\omega^{2}+4\gamma^{2}}}\sin\left(\frac{1}{2}t\sqrt{\omega^{2}+4\gamma^{2}}\right)\right],$
(23)

and

$\displaystyle\psi_{2}(x,t)=\Xi_{2}(t)\int_{-\infty}^{\infty}dx^{\prime}{\cal{G}}(x-x^{\prime},t)\varphi_{1}(x^{\prime})\;.$
(24)

The function $\Xi_{2}(t)$ is written as follows:

$\displaystyle\Xi_{2}(t)=-\frac{2i\gamma}{\sqrt{\omega^{2}+4\gamma^{2}}}e^{-\frac{i}{2}\omega
t}\sin\left(\frac{1}{2}t\sqrt{\omega^{2}+4\gamma^{2}}\right),$ (25)

and ${\cal{G}}(x,t)$ is the quantum free particle propagator, i.e.,
${\cal{G}}(x,t)=e^{-\frac{x^{2}}{2it}}/\sqrt{2\pi it}$.

In addition to the analytical result, given by the Eqs. (22) and (24), we also
obtain the numerical solutions of the Eqs. (15) and (16). For the numerical
approach, we consider the finite difference method [65]. We consider a grid
defined by $[0,X]\times[0,T]$, with boundary conditions equal to
$\psi_{1,2}(0,t)=\psi_{1,2}(X,t)=0$. The time is discretized by $t_{j}=j\Delta
t$, where $j=1,2,...,N_{t}$, with time step equal to $\Delta t=T/N_{t}$; and
the each space coordinate is given by $x_{i}=i\Delta x$, where
$i=1,2,...,N_{x}$, with space step equal to $\Delta x=X/N_{x}$. To avoid
numerical boundary problems, the origin of our space coordinate is in $X/2$.
From these considerations, the discretization of Eqs. (15) and (16) are given
by

$\psi^{i,j+1}_{1}=\psi^{i,j}_{1}+i\xi(\psi^{i+1,j}_{1}-2\psi^{i,j}_{1}+\psi^{i-1,j}_{1})-i\beta(V^{j}_{1}\psi^{i,j}_{2}+V^{j+1}_{1}\psi^{i,j+1}_{2}),$
(26)

and

$\psi^{i,j+1}_{2}=\psi^{i,j}_{2}+i\xi(\psi^{i+1,j}_{2}-2\psi^{i,j}_{2}+\psi^{i-1,j}_{2})-i\beta(V^{j}_{2}\psi^{i,j}_{1}+V^{j+1}_{2}\psi^{i,j+1}_{1}),$
(27)

where $\xi\equiv\Delta t/2\Delta x^{2}$, $\beta\equiv\gamma\Delta t/2$,
$V_{1}=e^{i\omega t}$ and $V_{2}=e^{-i\omega t}$. For the stability
conditions, it is required that $\xi\leq 1/2$ and the $\beta$ order less than
$\xi$ order [65].

Considering
$\psi_{1}(x,0)=e^{-\frac{x^{2}}{2\sigma^{2}}}/(2\pi\sigma^{2})^{1/4}$ and
$\psi_{2}(x,0)=0$ as the initial condition, the results for $|\psi_{1}|^{2}$
and $|\psi_{2}|^{2}$ are displayed in Figs. 2(a) and (b), respectively. The
parameters considered in this simulation are: $\xi=0.0016$, $\gamma=0.5$,
$\Delta x=0.25$, $\Delta t=0.0002$, and $\sigma=0.4$. As observed in the
results without kinetic terms, the system starts mostly in $\psi_{1}$.
However, a transition occurs to $\psi_{2}$ state due to the external field.
This effect is present in the presence of kinetic terms. Numerically, we
observed that
$\int_{-\infty}^{\infty}dx(|\psi_{1}(x,t)|^{2}+|\psi_{2}(x,t)|^{2})=1$. It is
worth mentioning that if we decrease $\Delta x$, the oscillations due to the
potential become smoother.

Figure 2: Probability distribution with a kinetic term of finding the system
in $\psi_{1}$ state, in the panel (a), and in $\psi_{2}$ state, in the panel
(b). The initial condition is given by
$\psi_{1}(x,0)=e^{-\frac{x^{2}}{2\sigma^{2}}}/(2\pi\sigma^{2})^{1/4}$ and
$\psi_{2}(x,0)=0$. We consider $\xi=0.0016$, $\gamma=0.5$, $\sigma=0.4$,
$\omega=2\pi$, $\Delta x=0.25$, $\Delta t=0.0002$, and $\sigma=0.4$.

The probability of finding the system in both states becomes approximately
equal after $t\geq 10$. The first state is mostly populated for a short time,
as observed in the result in Fig. 3. This result shows that the package
centered in origin is spread in the space in the first state and starts
transit to the second state in a sinusoidal form.

Figure 3: Probability distribution at $t=1.5$, for the states $\psi_{1}$ (red
line) and $\psi_{2}$ (blue line). We consider $\xi=0.0012$, $\gamma=0.5$,
$\sigma=0.4$, $\omega=2\pi$, $\Delta x=0.2$, $\Delta t=0.0001$, and
$\sigma=0.4$.

The mean square displacement is a measure of the spreading of the system,
represented by the wave package. It is widely applied in diffusion processes
to characterize the type of diffusion, usual or anomalous. For the usual
diffusion, we have a linear time dependence for the mean square displacement,
i.e., $\langle(\Delta x)^{2}\rangle\sim t$, which is related to the Markovian
processes. For the anomalous case, we have that $\langle(\Delta
x)^{2}\rangle\sim t^{S_{d}}$, where $S_{d}>1$ and $S_{d}<1$ are related to the
super-diffusive and sub-diffusive cases [3], respectively. In quantum
mechanics, we can also use this quantity to understand the spreading of the
probability density, i.e., $|\psi_{1,2}|^{2}$, in time. The normal case
corresponds to the free particle for the standard Schrödinger equation, where
$\langle(\Delta x)^{2}\rangle\sim t^{S}$ with $S=2$. The anomalous cases are
those that have different behaviors for the mean square displacement. Note
that these results are in agreement with the analytical results obtained for
Eq. (15) and (16), which results in Gaussian distributions for both wave
functions.

Considering the Gaussian package as the initial condition, the mean square
displacement for the free particle is shown in Fig. 4(a) by the black points,
which follows $\sim t^{S_{1}}$, with $S_{1}=2.02$. This result is obtained by
taking $\gamma=0$ in the numerical simulations. The effect of the potential is
displayed in Fig. 4(a) by the red points, which follows $\sim t^{S_{2}}$ with
$S_{2}=2.07$, for $|\psi_{1}|^{2}$. Due to the external potential, after a
certain time, the probability of finding the system transfer from the first
level to the second one, as shown in Fig. 4(b). The distribution for
$|\psi_{2}|^{2}$ increase as $\sim 2$. The slopes found by the numerical
simulations are in agreement with our analytical expressions, which indicate
$\sim t^{2}$ for both cases, free-particle and two-level system.

Figure 4: Mean square displacement for the Gaussian package. Panel (a) is for
$\psi_{1}$ state, and panel (b) is for $\psi_{2}$ standard case. The red
points are for standard two-level equations. The black points are for the free
particle. The slopes are $S_{1}=2.02$, $S_{2}=2.07$. We consider $\xi=0.0012$,
$\gamma=0.5$, $\sigma=0.4$, $\omega=2\pi$, $\Delta x=0.2$, $\Delta t=0.0001$,
and $\sigma=0.4$.

##  III Fractional Schrödinger Equations

Now, we analyze the previous scenarios within the fractional extension in time
of the Schrödinger equation. For the first case, i.e., the Hamiltonian given
by Eq. (4), we have that

$\displaystyle i^{\alpha}\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi_{1}(t)$ $\displaystyle=$ $\displaystyle
E_{1}\psi_{1}(t)+\gamma e^{i\omega t}\psi_{2}(t),$ (28)

and

$\displaystyle i^{\alpha}\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi_{2}(t)$ $\displaystyle=$ $\displaystyle
E_{2}\psi_{2}(t)+\gamma e^{-i\omega t}\psi_{1}(t),$ (29)

where $\hbar_{\alpha}=1$, without loss of generality. The relations
represented by Eqs. (28) and (29) extend Rabi’s model. For this set of
equations, obtaining an analytical solution for $\omega=0$, i.e., a static
field, is possible. To obtain the solutions for this case, we can use the
Laplace transform (${\cal{L}}\\{\psi(t);s\\}=\hat{\psi}(s)$ and
${\cal{L}}^{-1}\\{\hat{\psi}(s);t\\}=\psi(t)$) to simplify Eq. (28) and (29)
for the static case, yielding

$\displaystyle\hat{\psi}_{1}(s)=\frac{i^{\alpha}s^{\alpha-1}\left(i^{\alpha}s^{\alpha}-E_{2}\right)}{\left(i^{\alpha}s^{\alpha}-E_{1}\right)\left(i^{\alpha}s^{\alpha}-E_{2}\right)-\gamma^{2}},$
(30)

and

$\displaystyle\hat{\psi}_{2}(s)=\frac{i^{\alpha}s^{\alpha-1}}{\left(i^{\alpha}s^{\alpha}-E_{1}\right)\left(i^{\alpha}s^{\alpha}-E_{2}\right)-\gamma^{2}},$
(31)

for the initial condition $\psi_{1}(0)=1$ and $\psi_{2}(0)=0$. After
performing the inverse of the Laplace transform, it is possible to show that

$\displaystyle\psi_{1}(t)=\frac{1}{2}\left(1-\frac{E_{1}-E_{2}}{\sqrt{(E_{1}-E_{2})^{2}+4\gamma^{2}}}\right)E_{\alpha}\left(\gamma_{-}t^{\alpha}/i^{\alpha}\right)+\frac{1}{2}\left(1+\frac{E_{1}-E_{2}}{\sqrt{(E_{1}-E_{2})^{2}+4\gamma^{2}}}\right)E_{\alpha}\left(\gamma_{+}t^{\alpha}/i^{\alpha}\right),$
(32)

and

$\displaystyle\psi_{2}(t)=\frac{\gamma}{\sqrt{(E_{1}-E_{2})^{2}+4\gamma^{2}}}\bigg{(}E_{\alpha}\left(\gamma_{+}t^{\alpha}/i^{\alpha}\right)-E_{\alpha}\left(\gamma_{-}t^{\alpha}/i^{\alpha}\right)\bigg{)}\;,$
(33)

where
$\gamma_{\pm}=\left(E_{1}+E_{2}\pm\sqrt{(E_{1}-E_{2})^{2}+4\gamma^{2}}\right)/2$
and $E_{\alpha}(x)$ is the Mittag-Leffler function,

$\displaystyle E_{\alpha}(x)=\sum_{n=0}^{\infty}\frac{x^{n}}{\Gamma(1+\alpha
n)}\;,$ (34)

which corresponds to an extension of the exponential function [3]. The
solutions found for $\psi_{1}(t)$ and $\psi_{2}(t)$, are determined in terms
of the Mittag-Leffler function, implying that the system has an unusual
oscillation process, i.e., different from the standard case. For the case
$\omega\neq 0$, the solution can also be found, and it is given by

$\displaystyle\psi_{1}(t)$ $\displaystyle=$ $\displaystyle
E_{\alpha}\bigg{[}\left(E_{1}/i^{\alpha}\right)t^{\alpha}\bigg{]}$ (35)
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{\gamma}{i^{\alpha}}\right)^{2n}\int_{0}^{t}dt_{n}\Lambda(t-t_{n})\int_{0}^{t_{n}}dt_{n-1}\Lambda(t_{n}-t_{n-1})\cdots\int_{0}^{t_{2}}dt_{t_{1}}\Lambda(t_{2}-t_{1})E_{\alpha}\bigg{[}\left(E_{1}/i^{\alpha}\right)t_{1}^{\alpha}\bigg{]},$

with

$\displaystyle\psi_{2}(t)=\frac{\gamma}{i^{\alpha}}\int_{0}^{t}dt^{\prime}t^{\prime\alpha-1}E_{\alpha,\alpha}\bigg{[}\left(E_{2}/i^{\alpha}\right)(t-t^{\prime})^{\alpha}\bigg{]}e^{i\omega
t^{\prime}}\psi_{1}(t^{\prime})\;,$ (36)

where

$\displaystyle\Lambda(t)=e^{i\omega
t}\int_{0}^{t}dt^{\prime}t^{\prime\alpha-1}e^{-i\omega
t^{\prime}}E_{\alpha,\alpha}\bigg{[}\left(E_{1}/i^{\alpha}\right)t^{\prime\alpha}\bigg{]}(t-t^{\prime})^{\alpha-1}E_{\alpha,\alpha}\bigg{[}\left(E_{2}/i^{\alpha}\right)(t-t^{\prime})^{\alpha}\bigg{]}\;,$
(37)

by considering $\psi_{1}(0)=1$ and $\psi_{2}(0)=0$. The solutions for this
case are found in terms of the generalized Mittag-Leffler function [3],

$\displaystyle
E_{\alpha,\beta}(x)=\sum_{n=0}^{\infty}\frac{x^{n}}{\Gamma(\beta+\alpha
n)}\;.$ (38)

Figure 5 displays the numerical solution of Eqs. (28) and (29). For the static
case, in Figs. 5(a) and 5(b), and for the non-static case, in Figs. 5(c) and
5(d). The results are in perfect agreement with the analytical solutions found
in Eqs. (32), (33), (35), and (36) (see the Appendix for details of the
numerical procedure). A direct consequence obtained by incorporating
fractional time derivative in the Schrödinger equation is the non-conservation
of the probability, i.e.,
$|\psi_{1}(\infty)|^{2}+|\psi_{2}(\infty)|^{2}=1/\alpha^{2}$. This result
agrees with the results presented in Refs. [48, 47].

Figure 5: Probability of finding the system in $\psi_{1}(t)$ state given by
the red line, in panels (a) and (c), and in the $\psi_{2}(t)$ by the blue
line, in panels (b) and (d). The panel (a) is for $\omega=0$ and (b) for
$\omega=1$. We consider $E_{1}=1$, $E_{2}=2$, $\gamma=1$.

Considering the kinetic terms, the time fractional Schrödinger equation can be
written in the form

$i^{\alpha}\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi_{1}(x,t)=-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\psi_{1}(x,t)+\gamma e^{i\omega t}\psi_{2}(x,t),$ (39)

and

$i^{\alpha}\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi_{2}(x,t)=-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\psi_{2}(x,t)+\gamma e^{-i\omega t}\psi_{1}(x,t).$ (40)

These equations can be approximated by the following discretization [66]:

$\displaystyle\psi^{i,j+1}_{1}$ $\displaystyle=$
$\displaystyle\psi^{i,j}_{1}-i^{-\alpha}\xi_{\alpha}(\psi^{i+1,j}_{1}-2\psi^{i,j}_{1}+\psi^{i-1,j}_{1})+i^{-\alpha}\beta_{\alpha}(V^{j}_{1}\psi^{i,j}_{2}+V^{j+1}_{1}\psi^{i,j+1}_{2})$
(41) $\displaystyle-$
$\displaystyle\sum_{k=1}^{j}[(k+1)^{(1-\alpha)}-k^{(1-\alpha)}][\psi^{i,j+1-k}_{1}-\psi^{i,j-k}_{1}],$

and

$\displaystyle\psi^{i,j+1}_{2}$ $\displaystyle=$
$\displaystyle\psi^{i,j}_{2}-i^{-\alpha}\xi_{\alpha}(\psi^{i+1,j}_{2}-2\psi^{i,j}_{2}+\psi^{i-1,j}_{2})+i^{-\alpha}\beta_{\alpha}(V^{j}_{2}\psi^{i,j}_{1}+V^{j+1}_{2}\psi^{i,j+1}_{1})$
(43) $\displaystyle-$
$\displaystyle\sum_{k=1}^{j}[(k+1)^{(1-\alpha)}-k^{(1-\alpha)}][\psi^{i,j+1-k}_{2}-\psi^{i,j-k}_{2}],$

where $\xi_{\alpha}\equiv\Gamma(2-\alpha)\Delta t^{\alpha}/2\Delta x^{2}$,
$\beta_{\alpha}\equiv\Gamma(2-\alpha)\gamma\Delta t^{\alpha}/2$,
$V_{1}=e^{i\omega t}$ and $V_{2}=e^{-i\omega t}$. The convergence condition is
$\Delta t^{\alpha}/\Delta x^{2}\leq(1-2^{-\alpha})/\Gamma(2-\alpha)$ [67].

Figures 6(a) and 6(b) show the numerical solution for $\psi_{1}(x,t)$ and
$\psi_{2}(x,t)$ with $\alpha=0.98$, by considering the initial conditions
$\psi_{1}(x,0)=e^{-x^{2}/(2\sigma^{2})}/(2\pi\sigma^{2})^{1/4}$ and
$\psi_{2}(x,0)=0$, where $\sigma=0.4$. Note that, for $\alpha$ slightly
different from the standard case the dynamics properties of probabilities
densities spreads have a significant change from the standard case. If we
consider $\alpha<0.98$ these changes will be pronounced. For example, the
results presented in Fig. 6(a) and 6(b) show that the time fractional operator
makes the probability spread slowly when compared to the standard case. Also,
the transition between the states occurs with a greater amplitude than in the
integer case. Another anomalous behavior is non-probability conservation. In
this case, the probability decays, and the imaginary part of the effective
potential operates as a dissipate term [69].

Figure 6: Probability distribution with a kinetic term of finding the system
in $\psi_{1}$ state, in the panel (a), and in $\psi_{2}$ state, in the panel
(b), in the time fractional approach. We consider $\alpha=0.98$,
$\xi_{\alpha}=0.0012$, $\beta_{\alpha}=10^{-4}$, $\gamma=0.5$, $\sigma=0.4$,
$\omega=2\pi$, $\Delta x=0.5$, $\Delta t=0.0005$.

For a fixed time, a comparison between the probability distribution at the
space in the case where $\alpha=1$ (dotted lines) and $\alpha=0.98$
(continuous lines) is shown in Fig. 7. For this time, the results show that
the package spread decreases the amplitude of $|\psi_{1}|^{2}$, and the shape
of $|\psi_{2}|^{2}$ is wider.

Figure 7: Probability distribution at $t=1.5$, for the states $\psi_{1}$ (red
line) and $|\psi_{2}|^{2}$ (blue line). The continuous line shows the behavior
for $\alpha=0.98$ and the dotted for $\alpha=1$. We consider $\alpha=0.98$,
$\xi_{\alpha}=0.0012$, $\beta_{\alpha}=10^{-4}$, $\gamma=0.5$, $\sigma=0.4$,
$\omega=2\pi$, $\Delta x=0.5$, $\Delta t=0.0005$.

The mean square displacement for the Gaussian package is exhibited in Fig. 8
for $\alpha=0.98$ by the blue, in panel 8(a) for $|\psi_{1}|^{2}$ and in panel
8(b) for $|\psi_{2}|^{2}$. The red points are for the standard case, and the
black line is for the free particle. The fractional time spread is similar to
the standard case for short times. However, after this initial time, the blue
points follow $\sim t^{S_{3}}$, with $S_{3}=1.87$, while the red ones $\sim
t^{2.07}$. The behavior for the fractional case in time shows that the package
spreads with less intensity than the standard case; the spread is more
centered. The effect of the oscillatory potential is observed in the spread
for the second state, as shown in Fig. 8(b), by the blue curve for
$|\psi_{2}|^{2}$. The second state populated in fractional presence in this
scenario differs from the standard case. The fractional operator in time makes
the probability for $\psi_{2}$ state oscillates like a sinusoidal function. As
the time increase, the $\psi_{2}$ becomes populated with more frequency.
Another difference in fractional cases is that the probability is non-
conservative, and the deviations go to zero. The imaginary part of the
effective potential operates like a dissipate term [69].

Figure 8: Mean square displacement for the Gaussian package. Panel (a) is for
$\psi_{1}$ state, and panel (b) is for $\psi_{2}$ state. The blue points are
for $\alpha=0.98$, the red for standard equations ($S_{2}=2.07$), and the
black for a free particle. The slope associated with $\alpha=0.98$ is
$S_{3}=1.87$. We consider $\alpha=0.98$, $\xi_{\alpha}=0.0012$,
$\beta_{\alpha}=10^{-4}$, $\gamma=0.5$, $\sigma=0.4$, $\omega=2\pi$, $\Delta
x=0.5$, $\Delta t=0.0005$.

Now, let us consider the Schrödinger equation with fractional differential
operators in space. This extension can be directly related to the works of
Laskin [26], which takes Lévy flights in the Feynman path integral approach
into account. Following an analogous scheme [22], it is possible to include
the fractional differential operator in space in such a way that the equations
become

$i\frac{\partial}{\partial
t}\psi_{1}(x,t)=-\frac{1}{2}\frac{\partial^{\mu}}{\partial|x|^{\mu}}\psi_{1}(x,t)+\gamma
e^{i\omega t}\psi_{2}(x,t),$ (44)

and

$i\frac{\partial}{\partial
t}\psi_{2}(x,t)=-\frac{1}{2}\frac{\partial^{\mu}}{\partial|x|^{\mu}}\psi_{2}(x,t)+\gamma
e^{-i\omega t}\psi_{1}(x,t).$ (45)

This extension for the set of Schrödinger equation essentially considers
$\partial_{x}^{2}(\cdots)\rightarrow\partial_{|x|}^{\mu}(\cdots)$ with
$1<\mu<2$, which corresponds to a Riesz-Weyl fractional operator. By applying
the Fourier transform in the previous set of equations and using the property
${\mathcal{F}}\left\\{\partial_{|x|}^{\mu}\psi_{1,2}(x,t);k\right\\}=-|k|^{\mu}\widetilde{\psi}_{1,2}(k,t)$,
we have that

$i\frac{\partial}{\partial
t}\widetilde{\psi}_{1}(k,t)=\frac{1}{2}|k|^{\mu}\widetilde{\psi}_{1}(k,t)+\gamma
e^{i\omega t}\widetilde{\psi}_{2}(k,t),$ (46)

and

$i\frac{\partial}{\partial
t}\widetilde{\psi}_{2}(k,t)=\frac{1}{2}|k|^{\mu}\widetilde{\psi}_{2}(k,t)+\gamma
e^{-i\omega t}\widetilde{\psi}_{1}(k,t)\;.$ (47)

By performing some calculations, it is possible to show that

$\displaystyle i\frac{\partial}{\partial
t}\widetilde{\psi}_{1}(k,t)=\frac{1}{2}|k|^{\mu}\widetilde{\psi}_{1}(k,t)-i\gamma^{2}\int_{0}^{t}dt^{\prime}e^{-\frac{1}{2}i|k|^{\mu}(t-t^{\prime})}e^{i\omega(t-t^{\prime})}\widetilde{\psi}_{1}(k,t^{\prime})\;,$
(48)

which can be solved by using the Laplace transform.

The wave functions for this case can be obtained and are written as

$\displaystyle\widetilde{\psi}_{2}(k,t)=-i\gamma\int_{0}^{t}dt^{\prime}e^{-\frac{1}{2}i|k|^{\mu}(t-t^{\prime})}e^{-i\omega
t^{\prime}}\widetilde{\psi}_{1}(k,t^{\prime})\;,$ (49)

and

$\displaystyle\widetilde{\psi}_{1}(k,t)=e^{-\frac{1}{2}i\left(|k|^{\mu}-\omega\right)t}\left[\cos\left(\frac{1}{2}t\sqrt{\omega^{2}+4\gamma^{2}}\right)-\frac{i\omega}{\sqrt{\omega^{2}+4\gamma^{2}}}\sin\left(\frac{1}{2}t\sqrt{\omega^{2}+4\gamma^{2}}\right)\right]\widetilde{\varphi}_{1}(k)\;,$
(50)

assuming the initial conditions
$\widetilde{\psi}_{1}(k,0)=\widetilde{\varphi}_{1}(k)$ and
$\widetilde{\psi}_{2}(k,0)=0$. The inverse Fourier transform of Eqs. (49) and
(50) results in

$\displaystyle\psi_{1}(x,t)=\Xi_{1}(t)\int_{-\infty}^{\infty}dx^{\prime}{\cal{G}}_{\mu}(x-x^{\prime},t)\varphi_{1}(x^{\prime}),$
(51)

and

$\displaystyle\psi_{2}(x,t)=\Xi_{2}(t)\int_{-\infty}^{\infty}dx^{\prime}{\cal{G}}_{\mu}(x-x^{\prime},t)\varphi_{1}(x^{\prime}),$
(52)

with

${\cal{G}}_{\mu}(x,t)=\frac{1}{|x|}{\mbox{H}}^{1,1}_{2,2}\left[\frac{2}{it}|x|^{\mu}\left|{}^{\left(1,1\right),\left(1,\frac{\mu}{2}\right)}_{\left(1,\mu\right),\left(1,\frac{\mu}{2}\right)}\right.\right],$
(53)

which resembles the form of the Lévy distribution found in anomalous diffusion
processes. In Eq. (53), we have the H Fox function [70], usually represented
[3] by

$\displaystyle\mbox{\large{H}}^{m,n}_{p,q}\bigg{[}z\bigg{|}\begin{array}[]{c}\left(a_{p},A_{p}\right)\\\
\left(b_{q},B_{q}\right)\\\ \end{array}\bigg{]}$ $\displaystyle=$
$\displaystyle\mbox{\large{H}}^{m,n}_{p,q}\bigg{[}z\bigg{|}\begin{array}[]{c}\left(a_{1},A_{1}\right)\cdots\left(a_{p},A_{p}\right)\\\
\left(b_{1},B_{1}\right)\cdots\left(b_{q},B_{q}\right)\\\
\end{array}\bigg{]}=\frac{1}{2\pi i}\int_{L}ds\chi(s)z^{s}$ (58)

where

$\chi(s)=\frac{\prod_{j=1}^{m}\Gamma\left(b_{j}-B_{j}s\right)\prod_{j=1}^{n}\Gamma\left(1-a_{j}+A_{j}s\right)}{\prod_{j=1}^{q}\Gamma\left(1-b_{j}+B_{j}s\right)\prod_{j=1}^{p}\Gamma\left(a_{j}-A_{j}s\right)}\;,$
(59)

which involves Mellin–Barnes integrals [3]. The asymptotic behavior of Eq.
(53) in the limit of $|x|\rightarrow\infty$ is given by
${\cal{G}}_{\mu}(x,t)\sim i\left[t/\left(2|x|^{1+\mu}\right)\right]$, which is
different from the usual one characterized by the Gaussian behavior. Note that
this result for the asymptotic limit can be obtained by using the approach
employed in Ref. [71]. It is essentially an integration over the Mellin -
Branes integral poles, which represents Eq. (53). This feature is directly
connected to the presence of spatial fractional differential operators in Eqs.
(47) and (49).

In addition to the analytical approach, it is possible to investigate the
dynamical behavior from the numerical point of view by using the following
discretization

$\psi^{i,j+1}_{1}=\psi^{i,j}_{1}+i\xi_{\mu}\sum_{k=0}^{i-1}[\psi^{i-k+1,j}_{1}-2\psi^{i-k,j}_{1}+\psi^{i-k-1,j}_{1}][(k+1)^{2-\mu}-k^{2-\mu}]-i\beta(V^{j}_{1}\psi^{i,j}_{2}+V^{j+1}_{1}\psi^{i,j+1}_{2}),$
(60)

and

$\psi^{i,j+1}_{2}=\psi^{i,j}_{2}+i\xi_{\mu}\sum_{k=0}^{i-1}[\psi^{i-k+1,j}_{2}-2\psi^{i-k,j}_{2}+\psi^{i-k-1,j}_{2}][(k+1)^{2-\mu}-k^{2-\mu}]-i\beta(V^{j}_{2}\psi^{i,j}_{1}+V^{j+1}_{2}\psi^{i,j+1}_{1}),$
(61)

where $\xi_{\mu}\equiv\Delta t/2\Gamma(3-\mu)\Delta x^{\mu}$,
$\beta\equiv\gamma\Delta t/2$, $V_{1}=e^{i\omega t}$ and $V_{2}=e^{-i\omega
t}$ [72].

Small changes in the order of the fractional space operator make significant
changes in the spread probability dynamics. This phenomenon is observed in
Fig. 9, which exhibits the spread of $|\psi_{1}|^{2}$ in the panel (a) and
$|\psi_{2}|^{2}$ in the panel (b). The Gaussian package is the initial
condition spread widely in the presence of a space fractional operator
compared with the standard case. Another notable characteristic is the
behavior of $|\psi_{2}|^{2}$. The probability of finding the system in
$\psi_{2}$ state is more centered in origin and has a higher probability when
compared with the previous cases. The $\psi_{2}$ state assumes the Gaussian
shape probability and, for long times, replicates the dynamics observed in
$\psi_{1}$.

Figure 9: Probability distribution with a fractional kinetic term of finding
the system in $\psi_{1}$ state, in the panel (a), and in $\psi_{2}$ state, in
the panel (b). We consider $\mu=1.95$, $\xi_{\mu}=10^{-5}$, $\beta=10^{-5}$,
$\gamma=0.5$, $\sigma=0.4$, $\omega=2\pi$, $\Delta x=1.0$, $\Delta
t=0.000125$.

The comparison of the probabilities in $t=1.5$ is shown in Fig. 10, where the
continuous and dotted lines represent the cases $\mu=1.95$ and $\mu=2.0$,
respectively. The result shows a sharper division in the Gaussian package
along with an enlargement in the package. The probability of the $\psi_{2}$ is
more centered in origin, indicating that this state will take over a Gaussian
behavior over time.

Figure 10: Probability distribution at $t=1.5$ for the states $\psi_{1}$ (red
line) and $\psi_{2}$ (blue line). The continuous line shows the behavior for
$\mu=1.95$ and the dotted for $\mu=2$. We consider $\mu=1.95$,
$\xi_{\mu}=10^{-5}$, $\beta=10^{-5}$, $\gamma=0.5$, $\sigma=0.4$,
$\omega=2\pi$, $\Delta x=0.25$, $\Delta t=0.000125$.

The mean square displacement for the fractional derivative in space is shown
in Fig. 11, with the green points in Fig. 11(a) for $\psi_{1}$ and by the
green line in Fig. 11(b) for $\psi_{2}$. The red points are for the fractional
time derivative ($\alpha=0.98$), the blue points for the fractional space
derivative ($\mu=1.95$), and the black points for the free-particle. The
$|\psi_{1}|^{2}$ spread as $\sim t^{S_{4}}$ with $S_{4}=2.61$. Compared with
other cases, the fractional space operator makes the probability package
spread more widely, i.e., if we consider $\mu=2$ in a determined time, the
package occupies a certain range of space; however, for the same situation
with $\mu=1.95$ the package occupies a larger range. As shown in Fig. 11(b),
the $|\psi_{2}|^{2}$ spread more intensely. Populating the $\psi_{2}$ state
follows the Gaussian shape, as noted in Fig. 9(b).

Figure 11: Mean square displacement for the Gaussian package. Panel (a) is for
$\psi_{1}$ state and panel (b) is for $\psi_{2}$ state. The green points are
for $\mu=1.95$, the red for the standard case ($S_{2}=2.07$), the blue for
$\alpha=0.98$ ($S_{3}=1.87$), and the black for free-particle. The slope
associated with the green curve is $S_{4}=2.61$. We consider $\mu=1.95$,
$\xi_{\mu}=10^{-4}$, $\beta=10^{-5}$, $\gamma=0.5$, $\sigma=0.4$,
$\omega=2\pi$, $\Delta x=0.66$, $\Delta t=0.000125$.

The last possible case to be analyzed is the Schrödinger equation with
fractional differential operators in space and time by taking into account a
time-dependent potential, i.e.,

$i^{\alpha}\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi_{1}(x,t)=-\frac{1}{2}\frac{\partial^{\mu}}{\partial|x|^{\mu}}\psi_{1}(x,t)+\gamma
e^{i\omega t}\psi_{2}(x,t),$ (62)

and

$i^{\alpha}\frac{\partial^{\alpha}}{\partial
t^{\alpha}}\psi_{2}(x,t)=-\frac{1}{2}\frac{\partial^{\mu}}{\partial|x|^{\mu}}\psi_{2}(x,t)+\gamma
e^{-i\omega t}\psi_{1}(x,t)\;,$ (63)

It is possible to find a solution for these equations, and it is given by

$\displaystyle\psi_{1}(x,t)$ $\displaystyle=$
$\displaystyle\psi_{1}^{(0)}(x,t)+\sum_{n=1}^{\infty}\left(\gamma/i^{\alpha}\right)^{2n}\int_{-\infty}^{\infty}dx_{n}\int_{0}^{t}dt_{n}{\Upsilon}(x-x_{n},t-t_{n})$
(64) $\displaystyle\times$
$\displaystyle\int_{-\infty}^{\infty}dx_{n-1}\int_{0}^{t_{n}}dt_{n-1}{\Upsilon}(x_{n}-x_{n-1},t_{n}-t_{n-1})\cdots\int_{-\infty}^{\infty}dx_{1}\int_{0}^{t_{2}}dt_{1}{\Upsilon}(x_{2}-x_{1},t_{2}-t_{1})\psi_{1}^{(0)}(x_{1},t_{1})$

with

$\displaystyle\psi_{2}(x,t)=\frac{\gamma}{i^{\alpha}}\int_{-\infty}^{\infty}dx^{\prime}\int_{0}^{t}dt^{\prime}t^{\prime\alpha-1}{\cal{G}}_{\alpha,\mu}^{(\alpha)}(x-x^{\prime},t-t^{\prime})e^{i\omega
t^{\prime}}\psi_{1}(x^{\prime},t^{\prime})\;,$ (65)

where
$\psi_{1}^{(0)}(x,t)=\int_{-\infty}^{\infty}dx^{\prime}\varphi(x^{\prime}){\cal{G}}_{\alpha,\mu}^{(1)}(x-x^{\prime},t)$,

$\displaystyle\Upsilon(x,t)=e^{i\omega
t}\int_{-\infty}^{\infty}dx^{\prime}\int_{0}^{t}dt^{\prime}t^{\prime\alpha-1}e^{-i\omega
t^{\prime}}{\cal{G}}_{\alpha,\mu}^{(\alpha)}(x^{\prime},t^{\prime}){\cal{G}}_{\alpha,\mu}^{(\alpha)}(x-x^{\prime},t-t^{\prime})\;,$
(66)

and

$\displaystyle{\cal{G}}_{\alpha,\mu}^{(\beta)}(x,t)=\frac{1}{|x|}{\mbox{H}}^{2,1}_{2,3}\left[-\frac{|x|^{\mu}}{t^{\alpha}/(2i^{\alpha})}\left|{}^{\left(1,1\right),\left(\beta,\alpha\right),\left(1,\frac{\mu}{2}\right)}_{\left(1,\mu\right),\left(1,1\right),\left(1,\frac{\mu}{2}\right)}\right.\right]\;.$
(67)

Note that Eq. (67) is essentially the Green function of this case and,
consequently, connected to the relaxation process of this system. It differs
from the previous case since it mixes different fractional operators in space
and time.

It is possible to write a combination of the previous discretizations schemes
and find the equations

$\displaystyle\psi^{i,j+1}_{1}$ $\displaystyle=$
$\displaystyle\psi^{i,j}_{1}-\sum_{k=1}^{j}[(k+1)^{(1-\alpha)}-k^{(1-\alpha)}][\psi^{i,j+1-k}_{1}-\psi^{i,j-k}_{1}]+i^{-\alpha}\beta_{\alpha,\mu}(V^{j}_{1}\psi^{i,j}_{2}+V^{j+1}_{1}\psi^{i,j+1}_{2})$
(68) $\displaystyle-$ $\displaystyle
i^{-\alpha}\xi_{\alpha,\mu}\sum_{k=0}^{i-1}[\psi^{i-k+1,j}_{1}-2\psi^{i-k,j}_{1}+\psi^{i-k-1,j}_{1}][(k+1)^{2-\mu}-k^{2-\mu}],$

and

$\displaystyle\psi^{i,j+1}_{2}$ $\displaystyle=$
$\displaystyle\psi^{i,j}_{2}-\sum_{k=1}^{j}[(k+1)^{(1-\alpha)}-k^{(1-\alpha)}][\psi^{i,j+1-k}_{2}-\psi^{i,j-k}_{2}]+i^{-\alpha}\beta_{\alpha,\mu}(V^{j}_{2}\psi^{i,j}_{1}+V^{j+1}_{2}\psi^{i,j+1}_{1})$
(69) $\displaystyle-$ $\displaystyle
i^{-\alpha}\xi_{\alpha,\mu}\sum_{k=0}^{i-1}[\psi^{i-k+1,j}_{2}-2\psi^{i-k,j}_{2}+\psi^{i-k-1,j}_{2}][(k+1)^{2-\mu}-k^{2-\mu}],$

where $\xi_{\alpha,\mu}=\Gamma(2-\alpha)\Delta t^{\alpha}/2\Gamma(3-\mu)\Delta
x^{\mu}$, $\beta_{\alpha,\mu}=\gamma\Delta t^{\alpha}\Gamma(2-\alpha)/2$,
$V_{1}=e^{i\omega t}$ and $V_{2}=e^{-i\omega t}$. Considering $\alpha=0.98$
and $\mu=1.95$ the results for probability distribution is shown in Fig. 12(a)
for $|\psi_{1}|^{2}$ and in 12(b) for $|\psi_{2}|^{2}$. The results of the
combination of both fractional derivatives show a combination of the two
previously discussed behavior. However, for this set of parameters, the
results resemble fractional time dependence than space one.

Figure 12: Probability distribution with a fractional kinetic term of finding
the system in $\psi_{1}$ state, in the panel (a), and in $\psi_{2}$ state, in
the panel (b) for time and space fractional dependence. We consider
$\alpha=0.98$, $\mu=1.95$, $\gamma=0.5$, $\sigma=0.4$, $\omega=2\pi$,
$\xi_{\alpha,\mu}=6\times 10^{-4}$, $\beta_{\alpha,\mu}=\times 10^{-4}$,
$\Delta x=0.66$, $\Delta t=5\times 10^{-4}$.

Figure 13 displays the comparison between the probability distribution for
fractional time and space order (in continuous lines) versus the standard
model. These results make it possible to verify the composition of both
fractional operators. To gain greater influence for $\mu$ or $\alpha$ in the
dynamics, it is necessary to decrease one of these values.

Figure 13: Probability distribution at $t=1.5$ for the states $\psi_{1}$ (red
line) and $\psi_{2}$ (blue line). The continuous line show the behavior for
$\alpha=0.98$, $\mu=1.95$ and the dotted for $\alpha=1$ and $\mu=2$. We
consider $\alpha=0.98$, $\mu=1.95$, $\gamma=0.5$, $\sigma=0.4$, $\omega=2\pi$,
$\xi_{\alpha,\mu}=6\times 10^{-4}$, $\beta_{\alpha,\mu}=\times 10^{-4}$,
$\Delta x=0.66$, $\Delta t=5\times 10^{-4}$.

The spreading of the Gaussian package for $\alpha<1$ and $\mu<2$ is more
centered than in the standard case, as we see in Fig. 14(a) – orange points.
Figure 14(b) displays the spread for $\psi_{2}$ state by the orange line. The
associated slope is equal to $S_{5}=1.76$, near $S_{3}=1.87$, obtained for the
case when only the time fractional operator was considered. Therefore, in the
presence of both fractional operators, the time derivative supplants the
effects of the space derivative. The orange and blue points match in the range
displayed in Fig. 14(a) after some $t$. However, the distribution for
$|\psi_{2}|^{2}$, corresponding to the orange line in the panel 14(b), follows
the same shape as in the case of the fractional time derivative.

Figure 14: Mean square displacement for the Gaussian package. Panel (a) is for
$\psi_{1}$ state, and panel (b) is for $\psi_{2}$. The orange points are for
$\alpha=0.98$ and $\mu=1.95$, the green points for $\mu=1.95$, the red for
$\alpha=0.98$, and the black for a free particle. The slope associated with
$\alpha=0.98$ and $\mu=1.95$ is $S_{5}=1.76$. We consider $\alpha=0.98$,
$\mu=1.95$, $\gamma=0.5$, $\sigma=0.4$, $\omega=2\pi$,
$\xi_{\alpha,\mu}=6\times 10^{-4}$, $\beta_{\alpha,\mu}=\times 10^{-4}$,
$\Delta x=0.66$, $\Delta t=5\times 10^{-4}$.

##  IV Conclusion

We analyzed the influence of fractional operators in the Schrödinger equation
when an oscillating time-dependent potential is considered to simulate an
oscillatory external field applied in the system. We started with a two-level
system, which was first analyzed by considering the static case $\omega=0$ and
after the time-dependent case $\omega\neq 0$. We obtained time analytical and
numerical solutions for the standard and the fractional cases. In particular,
we verified that the solutions had an oscillating behavior for a long time.
Afterward, we incorporated the kinetic term in the Hamiltonian to allow the
spreading of the system. We also considered one state populated as an initial
condition while the other remained empty. We also analyzed this scenario from
the analytical and numerical points of view for the standard and the
fractional cases. For the fractional cases, we first consider the effect of
the fractional time derivatives, and after analyzing the spatial fractional
derivatives, which preserve the probability of the system. One of them is the
non-conservation of the probability of the system. We analyzed the behavior of
the mean square displacement (or deviation) for these cases and compared it
with the free particle case. The results showed that the fractional
differential operators lead to different behavior for spreading the system
when compared with the standard case. For the fractional derivative in space,
we have a faster spreading of the initial condition. On the other hand, we see
a slower spreading of the wave package when fractional derivatives in time are
incorporated in the Schrodinger equation. This feature is also present in the
diffusion context when fractional differential operators are considered,
evidencing that these operators strongly influence the random process
connected to these phenomena. The mean square displacement also evidenced this
point and the influence on the uncertain relations, as observed by Laskin
[25].

## Acknowledgements

The authors thank the financial support from the Brazilian Federal Agencies
(CNPq), the São Paulo Research Foundation (FAPESP, Brazil), CAPES, Fundação
Araucária. The authors thank the 105 Group Science (www.105groupscience.com).
E.K.L. acknowledges the support of the CNPq (Grant No. 301715/2022-0).

## Appendix I

A numerical method to solve initial-problem based on Caputo definition is a
generalization of the classical Adams-Bashforth-Moulton. This method was
proposed by Diethelm, Ford and Freed [73], and is defined by the follows
equations:

$y_{h}(t_{n+1})=\sum_{k=0}^{\left\lceil{\alpha}\right\rceil}\frac{t^{k}_{n+1}}{k!}y_{0}^{(k)}+\frac{h^{\alpha}}{\Gamma(\alpha+2)}f(t_{n+1},y^{P}_{h}(t_{n+1}))+\frac{h^{\alpha}}{\Gamma(\alpha+2)}\sum_{j=0}^{n}a_{j,n+1}f(t_{j},y_{h}(t_{j})),$
(70)

where

$y^{P}_{h}(t_{n+1})=\sum_{k=0}^{\left\lceil{\alpha}\right\rceil}\frac{t^{k}_{n+1}}{k!}y_{0}^{(k)}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{n}b_{j,n+1}f(t_{j},y_{h}(t_{j})).$
(71)

The coefficients are defined by

$a_{j,n+1}=\left\\{\begin{array}[]{lll}n^{\alpha+1}-(n-\alpha)(n+1)^{\alpha},\
\ {\rm if}\ j=0,\\\ (n-j+2)^{\alpha+1}+(n-j)^{\alpha+1}-2(n-j+1)^{\alpha+1},\
\ {\rm if}\ 1\leq j\leq n\\\ 1,\ \ {\rm if}j=n+1,\end{array}\right.$ (72)

and

$b_{j,n+1}=\frac{h^{\alpha}}{\alpha}((n+1-j)^{\alpha}-(n-j)^{\alpha}),$ (73)

where $j=1,2,...,n$ and $n$ is associated with discrete time window, $T$,
which is discrete in $t_{n}=nh$, with $n=0,1,...,N$, where $T=Nh$.

## References

  * [1] Gabrick EC, Sayari E, de Castro ASM, Trobia J, Batista AM, Lenzi EK. Fractional Schrödinger equation and time dependent potentials. Communications in Nonlinear Science and Numerical Simulation (2023), doi: https://doi.org/10.1016/j.cnsns.2023.107275
  * [2] Guo B, Pu X, Huang F. Fractional Partial Differential Equations and Their Numerical Solutions. World Scientific; 2015.
  * [3] Evangelista LR, Lenzi EK. Fractional Diffusion Equations and Anomalous Diffusion. Cambridge: Cambridge University Press; 2018.
  * [4] Solís-Pérez JE, Gómez-Aguilar JF, Atangana A. Novel numerical method for solving variable-order fractional differential equations with power, exponential and Mittag-Leffler laws. Chaos, Solitons and Fractals 2018;114:175-185.
  * [5] Herrmann R. Fractional calculus: an introduction for physicists. World Scientific; 2014.
  * [6] Templos-Hernández DJ, Quezada-Téllez LA, Gonzáles-Hernández BM, Rojas-Vite G, Pineda-Sánchez JE, Fernández-Anaya G, Rodriguez-Torres EE. A fractional-order approach to cardiac rhythm analysis. Chaos, Solitons and Fractals 2021;147:110942.
  * [7] Ciuchi F, Mazzulla A, Scaramuzza N, Lenzi EK, Evangelista LR. Fractional Diffusion Equation and the Electrical Impedance: Experimental Evidence in Liquid-Crystalline Cells. Journal Physics Chemistry C 2012;114:8773-8777.
  * [8] Bisquert J. Interpretation of a fractional diffusion equation with nonconserved probability density in terms of experimental systems with trapping or recombination. Physical Review E 2005;72:011109.
  * [9] Somer A, Novatski A, Serbena FC, Lenzi EK. Fractional GCEs behaviors merged: Prediction to the photoacoustic signal obtained with subdiffusive and superdiffusive operators. Journal of Applied Physics 2020;128:075107.
  * [10] Ali MS, Narayanan G, Shekher V, Alsaedi A, Ahmad B. Global Mittag-Leffler stability analysis of impulsive fractional-order complex-valued BAM neural networks with time varying delays. Communications in Nonlinear Science and Numerical Simulation 2020;83:105088.
  * [11] Pandey V, Holm S. A fractional calculus approach to the propagation of waves in an unconsolidated granular medium. The Journal of the Acoustical Society of America 2015;138(3):1766-1766.
  * [12] Bagley RL, Torvik PJ. A theoretical basis for the application of fractional calculus to viscoelasticity. Journal of Rheology 1983;27(3):201-210.
  * [13] Rosseto MP, Evangelista LR, Lenzi EK, Zola RS, Ribeiro de Almeida RR. Frequency-Dependent Dielectric Permittivity in Poisson–Nernst–Planck Model. The Journal of Physical Chemistry B 2022;126(34):6446–6453.
  * [14] Scarfone AM, Barbero G, Evangelista LR, Lenzi EK. Anomalous Diffusion and Surface Effects on the Electric Response of Electrolytic Cells. Physchem 2022;2(2):163-178.
  * [15] Lenzi EK, Guilherme LMS, da Silva BVHV, Koltun APS, Evangelista LR, Zola RS. Anomalous diffusion and electrical impedance response: Fractional operators with singular and non-singular kernels. Communications in Nonlinear Science and Numerical Simulation 2021;102:105907.
  * [16] Chen W, Hu S, Cai W. A causal fractional derivative model for acoustic wave propagation in lossy media. Archive of Applied Mechanics 2016;86(3):529-539.
  * [17] Cai W, Chen W, Fang J, Holm Sl. A survey on fractional derivative modeling of power-law frequency-dependent viscous dissipative and scattering attenuation in acoustic wave propagation. Applied Mechanics Reviews 2018;70(3):030802.
  * [18] Jiang Y, Qi H, Xu H, Jiang X. Transient electroosmotic slip flow of fractional Oldroyd-B fluids. Microfluidics and Nanofluidics 2017;21(1):1-10.
  * [19] Chang A, Sun HG, Zhang Y, Zheng C, Min F. Spatial fractional Darcy’s law to quantify fluid flow in natural reservoirs. Physica A: Statistical Mechanics and its Applications 2019;519:119-126.
  * [20] Chang A, Sun HG, Zheng C, Lu B, Lu C, Ma R, Zhang Y. A time fractional convection–diffusion equation to model gas transport through heterogeneous soil and gas reservoirs. Physica A: Statistical Mechanics and its Applications 2018;502:356-369.
  * [21] Pandey V, Holm S. Connecting the grain-shearing mechanism of wave propagation in marine sediments to fractional order wave equations. The Journal of the Acoustical Society of America 2016;140(6):4225-4236.
  * [22] Wang S, Xu M. Generalized fractional Schrödinger equation with space-time fractional derivatives. Journal of Mathematical Physics 2007;48:043502.
  * [23] Heydari MH, Razzaghi M, Baleanu D. A numerical method based on the piecewise Jacobi functions for distributed-order fractional Schrödinger equation. Communications in Nonlinear Science and Numerical Simulation 2023;116:106873.
  * [24] Laskin N. Fractional Quantum Mechanics. World Scientific Publishing Company; 2018.
  * [25] Laskin N. Fractional Schrödinger equation. Physical Review E 2002;66(5):056108.
  * [26] Laskin N. Fractals and quantum mechanics. Chaos: An Interdisciplinary Journal of Nonlinear Science 2000;10:780-790.
  * [27] Sandev T, Petreska I, Lenzi EK. Time-dependent Schrödinger-like equation with nonlocal term. Journal of Mathematical Physics 2014;55(9):092105.
  * [28] Lenzi EK, de Oliveira BF, da Silva LR, Evangelista LR. Solutions for a Schrödinger equation with a nonlocal term. Journal of Mathematical Physics 2008;49(3):032108.
  * [29] Sandev T, Petreska I, Lenzi EK. Generalized time-dependent Schrödinger equation in two dimensions under constraints. Journal of Mathematical Physics 2018;59(1):012104.
  * [30] Petreska I, de Castro ASM, Sandev T, Lenzi EK. The time-dependent Schrödinger equation in three dimensions under geometric constraints. Journal of Mathematical Physics 2019;60(3):032101.
  * [31] Sandev T, Petreska I, Lenzi EK. Constrained quantum motion in $\delta$-potential and application of a generalized integral operator. Computers & Mathematics with Applications 2019;78(5):1695-1704.
  * [32] Capelas de Oliveira E, Vaz Jr J. Tunneling in fractional quantum mechanics. Journal of Physics A: Mathematical and Theoretical 2011;44:185303.
  * [33] Guo X, Xu M. Some physical applications of fractional Schrödinger equation. Journal of Mathematical Physics 2006;47:082104.
  * [34] Dong J. Fractional Green’s Function for the Time-Dependent Scattering Problem in the Space-Time-Fractional Quantum Mechanics. International Journal of Theoretical Physics 2014;53:4065–4078.
  * [35] Naber M. Time fractional Schrödinger equation. Journal of Mathematical Physics 2004;45(8):3339-3352.
  * [36] Feynman RP, Hibbs AR. Quantum Mechanics and Path Integrals. New York:McGraw-Hill; 1965.
  * [37] Iomin A. Fractional-time Schrödinger equation: Fractional dynamics on a comb. Chaos, Solitons and Fractals 2011;44:348-352.
  * [38] Lenzi EK, Evangelista LR, Ribeiro HV, Magin RL. Schrödinger Equation with Geometric Constraints and Position-Dependent Mass: Linked Fractional Calculus Models. Quantum Reports 2022;4(3):296-308.
  * [39] Okposo NI, Veeresha A, Okposo EN. Solutions for time-fractional coupled nonlinear Schrödinger equations arising in optical solitons. Chinese Journal of Physics 2022;77:965-984.
  * [40] Achar BN, Narahari, Yale BT, Hanneken JW. Time Fractional Schrödinger Equation Revisited. Advances in Mathematical Physics 2013;2013:290216.
  * [41] Esen A, Sulaiman TA, Bulut H, Baskonus HM. Optical solitons to the space-time fractional (1+1)-dimensional coupled nonlinear Schrödinger equation. Optik 2018;167:150-156.
  * [42] Liaqat MI, Akgül A. A novel approach for solving linear and nonlinear time-fractional Schrödinger equations. Chaos, Solitons and Fractals 2022;162:112487.
  * [43] Hilfer R. Applications of Fractional Calculus in Physics. World Scientific; 2000.
  * [44] Heydari MH, Atangana A. A cardinal approach for nonlinear variable-order time fractional Schrödinger equation defined by Atangana–Baleanu–Caputo derivative. Chaos, Solitons and Fractals 2009;128:339-348.
  * [45] El-Nabulsi RA, Anukool W. A family of nonlinear Schrödinger equations and their solitons solutions. Chaos, Solitons and Fractals 2023;166:112907.
  * [46] Ain QT, He, J-H, Anjum N, Ali M. The fractional complex transform: a novel approach to the time-fractional Schrd̈inger equation. Fractals 2020;28(7):2050141.
  * [47] Zu C, Yu X. Time fractional Schrödinger equation with a limit based fractional derivative. Chaos, Solitons and Fractals 2022;157:111941.
  * [48] Lu L, Yu X. Time fractional evolution of the two-level system interacting with light field. Laser Physics Letters 2017;14(11):115202.
  * [49] Sakurai JJ, Napolitano J. Modern Quantum Mechanics. Cambridge:Cambridge University Press; 2017.
  * [50] Cohen-Tannoudji C, Diu B, Laloe F. Quantum Mechanics. Wiley-Interscience; 2006.
  * [51] Ruyten WM. Magnetic and optical resonance of two-level quantum systems in modulated fields. I. Bloch equation approach. Physical Review A 1990;42(7):4226-4245.
  * [52] Angelo RM, Wreszinski WF. Two-level quantum dynamics, integrability, and unitary NOT gates. Physical Review A 2005;72:034105.
  * [53] Cius D, Menon Jr L, dos Santos MAF, de Castro ASM, Andrade FM. Unitary evolution for a two-level quantum system in fractional-time scenario. Physical Review E 2022;106:054126.
  * [54] Itano WM, Bergquist JC, Bollinger JJ, Gilligian JM, Heinzen DJ, Moore FL. Quantum projection noise: Population fluctuations in two-level systems. Physical Review A 1993;47(5):3554-3570.
  * [55] Kibs OV, Slepyan GYa, Maksimenko SA, Hoffmann A. Matter Coupling to Strong Electromagnetic Fields in Two-Level Quantum Systems with Broken Inversion Symmetry. Physical Review Letters 2009;102(2):023601.
  * [56] Rabi II. Space quantization in a gyrating magnetic field. Physical Review 1937;51:652-654.
  * [57] Evangelista LR, Lenzi EK. An Introduction to Anomalous Diffusion and Relaxation. Springer Nature; 2023.
  * [58] Bayin SS. Definition of the Riesz derivative and its application to space fractional quantum mechanics. Journal of Mathematical Physics 2016;57:123501.
  * [59] Viñales AD, Despósito MA. Anomalous diffusion induced by a Mittag-Leffler correlated noise. Physical Review E 2007;75:042102.
  * [60] Despósito MA, Viñales AD. Memory effects in the asymptotic diffusive behavior of a classical oscillator described by a generalized Langevin equation. Physical Review E 2008;77:031123.
  * [61] Fa KS. Anomalous diffusion in a generalized Langevin equation. Journal of Mathematical Physics 2009;50:083301.
  * [62] Figueiredo Camargo R, Capelas de Oliveira E, Vaz Jr J. On anomalous diffusion and the fractional generalized Langevin equation for a harmonic oscillator. Journal of Mathematical Physics 2009;50:123518.
  * [63] Viñales AD, Wang KG, Despósito AM. Anomalous diffusive behavior of a harmonic oscillator driven by a Mittag-Leffler noise. Physical Review E 2009;80:011101.
  * [64] Figueiredo Camargo R, Chiacchio AO, Charnet R, Capelas de Oliveira E. Solution of the fractional Langevin equation and the Mittag–Leffler functions. Journal of Mathematical Physics 2009;50:063507.
  * [65] Crank J. The mathematics of Diffusion. Oxford University Press; 1975.
  * [66] Murio DA. Implicit finite difference approximation for time fractional diffusion equations. Computers and Mathematics with Applications 2008;56:1138-1145.
  * [67] Liu F, Shen S, Turner IW. Analysis of a Discrete non-Markovian Random Walk Approximation for the Time Fractional Diffusion Equation. The ANZIAM Journal 2005;46:C488-C504.
  * [68] Rydin YL, Mattsson K, Werpes J, Sjöqvist E. High-order ﬁnite difference method for the Schrödinger equation on deforming domains. Journal of Computational Physics 2021;443:110530.
  * [69] Bayin SS. Time fractional Schrödinger equation: Fox’s H-functions and the effective potential. Journal of Mathematical Physics 2013;54:012103.
  * [70] Mathai AM, Saxena RK, Haubold HJ. The H-Function. New York:Springer-Verlag New York; 2010.
  * [71] Saxena RK, Mathai AM, Haubold HJ. Fractional reaction-diffusion equations. Astrophysics and Space Science 2006;305:289-296.
  * [72] Shen S, Liu F. Error analysis of an explicit finite difference approximation for the space fractional diffusion equation with insulated ends. The ANZIAM Journal 2005;46:C871-C887.
  * [73] Diethelm K, Ford NJ, Freed AD, Luchko Y. Algorithms for the fractional calculus: A selection of numerical methods. Computer Methods in Applied Mechanics and Engineering 2005;194:743-773.