# Constraints on Late Time Violations of the Equivalence Principle in the Dark
Sector

Cameron C. Thomas and Carsten van de Bruck  Consortium for Fundamental
Physics, School of Mathematics and Statistics, University of Sheffield,
Hounsfield Road, Sheffield S3 7RH, United Kingdom

###### Abstract

If dark energy is dynamical due to the evolution of a scalar field, then in
general it is expected that the scalar is coupled to matter. While couplings
to the standard model particles are highly constrained by local experiments,
bounds on couplings to dark matter (DM) are only obtained from cosmological
observations and they are consequently weaker. It has recently be pointed out
that the coupling itself can become non-zero only at the time of dark energy
domination, due to the evolution of dark energy itself, leading to a violation
of the equivalence principle (EP) in the dark sector at late times. In this
paper we study a specific model and show that such late-time violations of the
EP in the DM sector are not strongly constrained by the evolution of the
cosmological background and by observables in the linear regime (e.g. from the
cosmic microwave background radiation). A study of perturbations in non-linear
regime is necessary to constrain late–time violations of the equivalence
principle much more strongly.

##  I Introduction

There are several reasons to study alternative theories to the cosmological
constant as a model for dark energy (DE). Firstly, if DE is due to a non-
vanishing cosmological constant, its value has to be very small to fit the
data. The expected theoretical value, however, is much larger. This problem
has not been solved, but there are attempts to address this problem (see
Weinberg (1989); Padmanabhan (2003); Martin (2012) and references therein).
Secondly, there are several tensions between data sets, providing tantalising
hints that the standard model of cosmology, the $\Lambda$-Cold-Dark-
Matter($\Lambda$CDM) model, may be in need of corrections (we refer to Abdalla
_et al._ (2022) for a recent overview over the tensions and Di Valentino  _et
al._ (2021) for an overview of attempts to solve the tensions). Among the
extensions of the $\Lambda$CDM model which remain the best motivated ones are
scalar-field models of DE, in which the accelerated expansion is driven by a
scalar field Wetterich (1988); Ratra and Peebles (1988); Caldwell  _et al._
(1998). It is expected that, in general, the scalar is coupled to at least one
species of matter, unless there is a symmetry which forbids such couplings
Carroll (1998). Such a coupling results in an additional force mediated
between the coupled species. Since the interaction between DE and ordinary
matter is strongly constrained, in some models only the coupling to cold dark
matter (CDM) is of cosmological significance Wetterich (1995); Amendola
(2000). It is this type of theory which we consider in this paper.

It has recently been suggested, that the potential energy of scalar fields
appearing in string theory cannot be arbitrarily flat Obied  _et al._ (2018);
Agrawal  _et al._ (2018), see Palti (2019) for an overview of the swampland
programme. If true, the accelerated expansion cannot be driven by a
cosmological constant (de Sitter space is not realised in string theory) and
the equation of state of dark energy is not constant and deviates potentially
significantly from the value expected in the $\Lambda$CDM model. Additionally,
a coupling of the scalar to some sectors in the theory are expected. Based on
these observations, several phenomenological models have been proposed
recently van de Bruck and Thomas (2019); Agrawal  _et al._ (2021); Olguin-
Trejo  _et al._ (2019); Brax  _et al._ (2020); Valeixo Bento  _et al._ (2020);
Brinkmann  _et al._ (2022). In this paper we study specific a model in which
the coupling function between the dark energy field and dark matter has a
minimum van de Bruck and Thomas (2019). As a result of the minimum, the
coupling switches on only at late times, at the beginning of the dark energy
dominated epoch. One of our main results of this paper is that the regime in
which linear perturbation theory is valid does not constrain the parameter of
the model greatly. In other words, late time violations of the equivalence
principle in the dark sector are not strongly constrained by studying the
background evolution or CMB anisotropies. Instead, to obtain stronger
constraints a study of the non-linear regime in considerable detail is needed.

The paper is organised as follows. In Section II we present the model and its
parameter. In Section III we describe our methodology, describe the data sets
used and present the constraints on the model. We conclude in Section IV.

##  II Interacting Dark Energy

The model we consider consists of the gravitational sector described by the
Einstein–Hilbert action without cosmological constant, a part which describes
the standard model (SM) particles and a part for DE described by a scalar
field $\phi$ with potential energy $V(\phi)$. Finally, the interaction between
DE and DM is described by a conformal coupling. The full action reads

$\displaystyle{\cal S}$ $\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-g}\left(\frac{M_{\rm Pl}^{2}}{2}{\cal
R}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-V(\phi)\right)$
(1) $\displaystyle+$ $\displaystyle\int d^{4}x\sqrt{-g}{\cal L}_{\rm
SM}(g,\Psi_{i})+\int d^{4}x\sqrt{-{\tilde{g}}}{\cal L}_{\rm
DM}({\tilde{g}},\chi)$

where $M_{\rm Pl}$ is the reduced Planck mass, ${\cal R}$ is the Ricci–scalar,
the SM fields are denoted by $\Psi_{i}$ and $\chi$ is the DM field (assuming
here for simplicity that dark matter consists of only one species). The
metrics $g$ and $\tilde{g}$ are related by a conformal transformation
$\tilde{g}_{\mu\nu}=C(\phi)g_{\mu\nu}$. Such theories have been discussed in
considerable length in the literature, but the new ingredient in this paper is
that the function $C(\phi)$ has a minimum. That the coupling functions in
string theory could possess a minimum due to non-perturbative effects was
suggested in the works by Damour and Polyakov Damour and Polyakov (1994) and
has been used in Brax  _et al._ (2010) to construct a dilaton-model of DE.
Combining these theoretical developments motivated us in van de Bruck and
Thomas (2019) to consider a specific model in which $C(\phi)$ has a minimum at
some value scalar field value $\phi_{*}$. To be concrete, in this paper we
consider the following function for $C(\phi)$:

$C(\phi)=\cosh{\left(\sqrt{\alpha}\,(\phi-\phi_{*})/M_{\rm Pl}\right),}$ (2)

where $\phi_{*}$ denotes the minimum of the function $C(\phi)$ and $\alpha$ is
a constant. In this paper we choose $\phi_{*}=1~{}M_{\rm Pl}$ without loss of
generality. In van de Bruck and Thomas (2019) it was pointed out that even if
the field starts away from the minimum in the very early universe, there are
attractor mechanism at work in the early universe which drives the field
towards the minimum quickly. Nevertheless, in our analysis below we allow the
field to start away from the minimum value at $\phi_{*}$ to find constraints
on the initial conditions of the DE field. In our analysis we choose an
exponential potential with

$V(\phi)=V_{0}e^{-\lambda\phi/M_{\rm Pl}},$ (3)

where $\lambda$ denotes the slope of the potential, which in string theory
according to Obied  _et al._ (2018) cannot be arbitrarily small and should be
${\cal O}(1)$. Finally, we assume in the following that the universe is
spatially flat.

Because of the coupling, there is an exchange of energy between DM and the DE
field. As a result, the evolution of the DM density and the modified
Klein–Gordon equation are given by

$\dot{\rho}_{\rm DM}+3H\rho_{\rm DM}=\beta M_{\rm Pl}^{-1}\dot{\phi}\rho_{\rm
DM}$

and

$\ddot{\phi}+3H\dot{\phi}+V_{\phi}=-\beta M_{\rm Pl}^{-1}\rho_{\rm DM}.$

In these equations we have defined

$\beta=\frac{M_{\rm Pl}}{2}\frac{{\rm d}\ln C}{{\rm d}\phi}.$

The effective gravitational constant between two DM particles is given by
Amendola (2004); van de Bruck and Thomas (2019)

$G_{\rm eff}=G_{N}\left(1+2\beta^{2}\right).$ (4)

The evolution of $G_{\rm eff}$ is shown in Figure 1 for various choices of the
parameter $\alpha$ and $\lambda$. In general, the additional force between DM
particles due to the scalar field only becomes significant at redshifts $z<1$,
when the DE field starts to evolve due to the influence of the potential,
thereby displacing it from the minimum of the coupling function. We emphasize
that the effective gravitational coupling (4) between DM particles depends on
$\alpha$ as well as on the scalar field.

To summarize, the parameter of the model we seek to constrain are the slope of
the potential $\lambda$, the parameter $\alpha$, which is related to strength
of the coupling between DM and DE and the initial field value $\phi_{\rm ini}$
deep inside the radiation dominated epoch.

Figure 1: The evolution of the effective gravitational constant, defined in
Equation 4, for models with a different value of $\alpha$ but same value of
$\lambda=0.1$ (top) and for models with a different value of $\lambda$ but
same value of $10^{-4}\alpha=1$ (bottom). Where applicable, the values of our
cosmological parameters are taken from the best fit values of a $\Lambda$CDM
cosmology based on $Planck$ TTTEEE+lowE, such as in Table 2 of Aghanim  _et
al._ (2020).

##  III Methodology, Data and Results

Table 1: Flat priors for the cosmological parameters sampled in our analysis. |
---|---
Parameter |  Prior
|
|
$\Omega_{b}h^{2}$ |  $[0.005,0.1]$
$\Omega_{cdm}h^{2}$ |  $[0.001,0.99]$
$100\theta_{s}$ |  $[0.5,10.0]$
$\ln\left(10^{10}A_{s}\right)$ |  $[2.7,4.0]$
$n_{s}$ |  $[0.9,1.1]$
$\tau_{reio}$ |  $[0.01,0.8]$
$\lambda$ |  $[0,5]$
$10^{-4}\alpha$ |  $[0,50]$
$\phi_{\rm ini}/M_{\rm Pl}$ |  $[-5,5]$
|

In our analysis, the IDE model is described by a set of nine parameters whose
priors are specified in Table 1 and where $h$ is the reduced Hubble constant
defined by $H_{0}=100h{\rm kms^{-1}Mpc^{-1}}$. These parameters are the
reduced baryon energy density $\Omega_{b}h^{2}$, the reduced CDM energy
density $\Omega_{cdm}h^{2}$, the ratio of the sound horizon to the angular
diameter distance at decoupling $\theta_{s}$, the scalar amplitude of the
primordial power spectrum $A_{s}$, the scalar spectral index $n_{s}$, the
reionization optical depth $\tau_{reio}$, the slope of the scalar field
potential $\lambda$, the conformal coupling parameter $\alpha$, and the
initial value of the scalar field $\phi_{\rm ini}$.

In order to numerically study the evolution of the background and cosmological
perturbations for the Interacting Dark Energy (IDE) model described above we
use a modified version of the CLASS code Lesgourgues (2011); Blas  _et al._
(2011) to calculate the background evolution and the evolution of
perturbations111https://github.com/lesgourg/class_public. For cosmological
parameter exploration we use the Markov Chain Monte Carlo sampling package
MontePython Brinckmann and Lesgourgues (2018); Audren  _et al._ (2013), in
conjunction with the data sets outlined
below222https://github.com/brinckmann/montepython_public. In addition to this,
the GetDist Lewis (2019) Python package333https://github.com/cmbant/getdist is
used to analyse the chains and produce the values and plots of the parameters
in Table 2 and Figure 2 .

We use the following combination of recent observational data sets in order to
analyse and constrain the IDE model:

  * •

Cosmic Microwave Background:

We use the full TTTEEE+lowE Cosmic Microwave Background (CMB) likelihood from
the latest Planck 2018 release Aghanim  _et al._ (2020). This includes
temperature (TT) and polarisation (EE) anisotropy data as well as cross-
correlation data between temperature and polarisation (TE) at high and low
multipoles.

  * •

Baryon Acoustic Oscillations:

We consider Baryon Acoustic Oscillations (BAO) measurements coming from BOSS
DR12 Alam  _et al._ (2017), 6dFGS Beutler  _et al._ (2011), and SDSS-MGS Ross
_et al._ (2015) for use in our analysis.

  * •

Type Ia Supernovae:

We use the Pantheon data catalog consisting of $1048$ points in the region
$z\in[0.01,2.3]$ of SNIa luminosity distance data as provided by Scolnic  _et
al._ (2018).

  * •

Redshift Space Distortions:

We employ the ’Gold 2018’ Redshift Space Distortions (RSD) data set
compilation consisting of 22 measurements as described in Sagredo  _et al._
(2018) and a likelihood code as detailed in Arjona  _et al._ (2020).

  * •

Cosmic Chronometers:

We use 31 measurements of $H(z)$ in the redshift range $z\in[0.07,1.965]$ as
detailed in Table 4 of Moresco  _et al._ (2016).

The results for our data analysis are shown in Table 2 where we compare
constraints on $\Lambda$CDM with the IDE model and the uncoupled model using
the full PL18+BAO+Pantheon+CC+RSD data set for all models, and also in Figure
2 where we plot marginalised posterior distributions for parameters in the IDE
model. We utilise the full PL18+BAO+Pantheon+CC+RSD data set in order to
obtain convergence for all the models, with all parameters achieving
$|R-1|<0.03$, where $R$ is the Gelman-Rubin statistic Gelman and Rubin (1992).
In addition to $\Lambda$CDM and the IDE model, we also consider the uncoupled
model (i.e. imposing $\alpha=0$) for comparison and to illustrate the impact
of including the coupling on the parameter constraints. For the uncoupled
model, owing to the form of the potential, we fix the initial value of the
scalar field, $\phi_{\rm ini}$, and do not allow it to vary in our analysis.

One of the first things to notice is the large values of the coupling
constant, $\alpha$, allowed by the data as can be seen in Figure 2. We are
only able to attain an upper limit for $\alpha$ of $10^{-4}\alpha<2.67$ at
$2\sigma$. As can be seen in the upper panel of Figure 1, the large values of
$\alpha$ supported by the data actually do not have as drastic an impact on
$G_{\rm eff}$ as one may be led to believe. We can see that a ten-fold
increase in $\alpha$ (from $10^{-4}\alpha=0.01$ to $10^{-4}\alpha=0.1$) leads
to an increase of $~{}5\%$ in $G_{\rm eff}$, but a ten-fold increase from
$10^{-4}\alpha=0.5$ to $10^{-4}\alpha=5$ leads to only a sub-percentage
increase in $G_{\rm eff}$. $G_{\rm eff}$ appears to be far more sensitive to
changes in the slope of the potential for the model, $\lambda$, as can be seen
in the lower panel of Figure 1. The reason for this is that $G_{\rm eff}$
depends on $\alpha$ and on the value of $\phi$. Larger values of $\alpha$
imply a steeper coupling function and as a consequence, the field is kept at
the minimum more effectively. For the most extreme cases, we find that the
field excursion is smaller than the Planck mass, i.e. ($|\phi-\phi_{*}|<M_{\rm
Pl}$).

We obtain an upper bound on the slope of the potential of $\lambda<0.192$ at
$2\sigma$ for the coupled IDE model and an upper bound of $\lambda<0.406$ at
$2\sigma$ for the uncoupled model. We can see that this IDE model does not
help to alleviate the swampland requirement of $|V^{\prime}/V|\geq
c\sim\mathcal{O}(1)$ and in fact exacerbates the requirement when compared to
the uncoupled model. We also note that the coupled model breaks the degeneracy
in the $(H_{0},\lambda)$-plane and the $(\sigma_{8},\lambda)$-plane, as it can
be seen in Figure 3.

From Figure 2, it is clear that there also exists a negative correlation
between the conformal coupling constant, $\alpha$, and slope of the scalar
field potential, $\lambda$, with lower values of $\alpha$ preferring higher
values of $\lambda$.

We justify the prior for the initial value of the scalar field, $\phi_{\rm
ini}$, to be centred around the minimum of the coupling function, $\phi_{*}$,
owing to the attractor mechanisms aforementioned. We allow $\phi_{\rm ini}$ to
vary as a cosmological parameter in our analysis and find that it is tightly
constrained about the minimum of the coupling function with a mean of
$\phi_{\rm ini}/M_{\rm Pl}=1.0004^{+0.0082}_{-0.0092}$, lending credence to
the attractor mechanism discussed in van de Bruck and Thomas (2019).

The IDE model is in strong agreement with $\Lambda$CDM for this data set
combination, as can be seen in Table 2, with no significant tensions observed
in the cosmological parameters between the two models. This is to be expected
as the model behaves like standard uncoupled quintessence (with an exponential
potential) up until $z\sim 1$, when the field switches on, at which point it
begins to deviate from an uncoupled model.

Figure 2: Marginalised posterior distributions for the IDE model using PL18+BAO+Pantheon+CC+RSD. Figure 3: Marginalised posterior distributions for the IDE and uncoupled quintessence ($\alpha=0$) models using PL18+BAO+Pantheon+CC+RSD. Table 2: Observational constraints at a $68\%$ confidence level on the independent and derived cosmological parameters for $\Lambda$CDM and the IDE model using PL18+BAO+Pantheon+CC+RSD. Where $68\%$ confidence limits could not be found, an upper $95\%$ confidence limit was used, i.e. for parameters $10^{-4}\alpha$ and $\lambda$. For comparison we also calculated the constraints for an uncoupled quintessence model (last column). The quantities in the second half of this table are the derived parameters of our analysis which are the redshift of reionization $z_{reio}$, the Helium fraction $Y_{He}$, the Hubble constant $H_{0}$, the absolute magnitude of SN1a as inferred from the data sets used $M$, and the present day mass fluctuation amplitude at $8h^{-1}{\rm Mpc}$ $\sigma_{8}$. |  |  |
---|---|---|---
Parameter |  $\Lambda$CDM |  IDE |  Uncoupled quintessence ($\alpha=0$)
|  |  |
$\Omega_{b}h^{2}$ |  $0.02247\pm 0.00013$ |  $0.02247\pm 0.00013$ |  $0.02248\pm 0.00013$
$\Omega_{cdm}h^{2}$ |  $0.11875\pm 0.00094$ |  $0.11864\pm 0.00098$ |  $0.11856\pm 0.00096$
$100\theta_{s}$ |  $1.04199\pm 0.00028$ |  $1.04202\pm 0.00028$ |  $1.04201\pm 0.00028$
$\ln\left(10^{10}A_{s}\right)$ |  $3.041\pm 0.016$ |  $3.041\pm 0.016$ |  $3.042\pm 0.016$
$n_{s}$ |  $0.9685\pm 0.0037$ |  $0.9687\pm 0.0039$ |  $0.9690\pm 0.0037$
$\tau_{reio}$ |  $0.0538\pm 0.0077$ |  $0.0540\pm 0.0080$ |  $0.0545\pm 0.0076$
$10^{-4}\alpha$ |  - |  $<2.67$ |  -
$\lambda$ |  - |  $<0.109$ |  $<0.406$
$\phi_{\rm ini}$ |  - |  $1.0004^{+0.0082}_{-0.0092}$ |  -
|  |  |
$z_{reio}$ |  $7.59\pm 0.78$ |  $7.60\pm 0.82$ |  $7.64\pm 0.77$
$Y_{He}$ |  $0.247880\pm 0.000057$ |  $0.247882\pm 0.000056$ |  $0.247885\pm 0.000057$
$H_{0}$ |  $67.96\pm 0.42$ |  $68.02\pm 0.44$ |  $67.67^{+0.56}_{-0.45}$
$M$ |  $-19.411\pm 0.012$ |  $-19.409\pm 0.012$ |  $-19.415\pm 0.013$
$\sigma_{8}$ |  $0.8055\pm 0.0068$ |  $0.8066\pm 0.0071$ |  $0.8019^{+0.0080}_{-0.0072}$
|  |  |

##  IV Summary and Conclusions

In this paper we have studied a specific IDE model, in which a fifth force
between DM particles switches on at the onset of DE domination. This is
achieved by considering a coupling function which possesses a minimum at a
certain field value $\phi_{*}$. The DE field evolves away from the minimum due
to the influence of the potential energy. In contrast to other interacting DE
models, the fifth force between DM particles becomes only important at the end
of matter domination, if the field does sit at the minimum of the coupling
function initially. We used several data sets to constrain $\phi_{*}$, the
coupling parameter $\alpha$ and the slope of the potential energy $\lambda$.

Our findings can be summarized as follows:

  * •

The initial value of the field, deep inside the radiation dominated epoch, is
constrained to be very near the minimum of the coupling function.

  * •

Our best fit value for the slope $\lambda$ of the potential is lower than
other IDE models, see e.g. van de Bruck and Mifsud (2018), although different
data sets have been used in this study.

  * •

The coupling parameter $\alpha$ is only weakly constrained
($10^{-4}\alpha\lesssim 2.6$).

  * •

Since the effective coupling depends also on the value of the DE field $\phi$,
the effective gravitational constant is determined by the dynamics of DE and
hence by both $\alpha$ and $\lambda$. The dependence on the effective
gravitational coupling on $\lambda$ and $\alpha$ is illustrated in 1.

Our study shows that equivalence principle violations in the dark sector at
the present epoch are much less constrained than previously thought. This is
because such new interactions may result from the dynamics of DE itself and
only become important rather late in the cosmic history. The data sets used in
this paper, which measure the evolution of the cosmological background and
cosmological perturbations at the linear level, do not constrain the fifth
force strongly and the model is, for a wide range of parameter,
indistinguishable from the $\Lambda$CDM model. We do expect the model to
differ from $\Lambda$CDM at much smaller length scales, for which linear
perturbation theory is no longer adequate. The fifth force is likely to leave
an imprint on the evolution of non-linear perturbations and affect, for
example, the matter power spectra on small scales. It would also be
interesting to study the tidal tail test Kesden and Kamionkowski (2006) in
this model, since the coupling switches on only at late times. We do expect
this test to constrain this model further. We leave these studies for future
work.

## Acknowledgements

We are grateful to Ed Copeland and Adam Moss for useful remarks during the
early stages of this work. CCT is supported by a studentship from the School
of Mathematics and Statistics at the University of Sheffield. CvdB is
supported (in part) by the Lancaster–Manchester–Sheffield Consortium for
Fundamental Physics under STFC grant: ST/T001038/1.

## References

  * Weinberg (1989) Steven Weinberg, “The Cosmological Constant Problem,” Rev. Mod. Phys. 61, 1–23 (1989).
  * Padmanabhan (2003) T. Padmanabhan, “Cosmological constant: The Weight of the vacuum,” Phys. Rept. 380, 235–320 (2003), arXiv:hep-th/0212290 .
  * Martin (2012) Jerome Martin, “Everything You Always Wanted To Know About The Cosmological Constant Problem (But Were Afraid To Ask),” Comptes Rendus Physique 13, 566–665 (2012), arXiv:1205.3365 [astro-ph.CO] .
  * Abdalla  _et al._ (2022) Elcio Abdalla _et al._ , “Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies,” JHEAp 34, 49–211 (2022), arXiv:2203.06142 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021) Eleonora Di Valentino, Olga Mena, Supriya Pan, Luca Visinelli, Weiqiang Yang, Alessandro Melchiorri, David F. Mota, Adam G. Riess, and Joseph Silk, “In the realm of the Hubble tension—a review of solutions,” Class. Quant. Grav. 38, 153001 (2021), arXiv:2103.01183 [astro-ph.CO] .
  * Wetterich (1988) C. Wetterich, “Cosmology and the Fate of Dilatation Symmetry,” Nucl. Phys. B 302, 668–696 (1988), arXiv:1711.03844 [hep-th] .
  * Ratra and Peebles (1988) Bharat Ratra and P. J. E. Peebles, “Cosmological Consequences of a Rolling Homogeneous Scalar Field,” Phys. Rev. D 37, 3406 (1988).
  * Caldwell  _et al._ (1998) R. R. Caldwell, Rahul Dave, and Paul J. Steinhardt, “Cosmological imprint of an energy component with general equation of state,” Phys. Rev. Lett. 80, 1582–1585 (1998), arXiv:astro-ph/9708069 .
  * Carroll (1998) Sean M. Carroll, “Quintessence and the rest of the world,” Phys. Rev. Lett. 81, 3067–3070 (1998), arXiv:astro-ph/9806099 .
  * Wetterich (1995) Christof Wetterich, “The Cosmon model for an asymptotically vanishing time dependent cosmological ’constant’,” Astron. Astrophys. 301, 321–328 (1995), arXiv:hep-th/9408025 .
  * Amendola (2000) Luca Amendola, “Coupled quintessence,” Phys. Rev. D 62, 043511 (2000), arXiv:astro-ph/9908023 .
  * Obied  _et al._ (2018) Georges Obied, Hirosi Ooguri, Lev Spodyneiko, and Cumrun Vafa, “De Sitter Space and the Swampland,” (2018), arXiv:1806.08362 [hep-th] .
  * Agrawal  _et al._ (2018) Prateek Agrawal, Georges Obied, Paul J. Steinhardt, and Cumrun Vafa, “On the Cosmological Implications of the String Swampland,” Phys. Lett. B 784, 271–276 (2018), arXiv:1806.09718 [hep-th] .
  * Palti (2019) Eran Palti, “The Swampland: Introduction and Review,” Fortsch. Phys. 67, 1900037 (2019), arXiv:1903.06239 [hep-th] .
  * van de Bruck and Thomas (2019) Carsten van de Bruck and Cameron C. Thomas, “Dark energy, the swampland, and the equivalence principle,” Physical Review D 100 (2019), 10.1103/physrevd.100.023515.
  * Agrawal  _et al._ (2021) Prateek Agrawal, Georges Obied, and Cumrun Vafa, “$H_{0}$ tension, swampland conjectures, and the epoch of fading dark matter,” Phys. Rev. D 103, 043523 (2021), arXiv:1906.08261 [astro-ph.CO] .
  * Olguin-Trejo  _et al._ (2019) Yessenia Olguin-Trejo, Susha L. Parameswaran, Gianmassimo Tasinato, and Ivonne Zavala, “Runaway Quintessence, Out of the Swampland,” JCAP 01, 031 (2019), arXiv:1810.08634 [hep-th] .
  * Brax  _et al._ (2020) Philippe Brax, Carsten van de Bruck, and Anne-Christine Davis, “Swampland and screened modified gravity,” Phys. Rev. D 101, 083514 (2020), arXiv:1911.09169 [hep-th] .
  * Valeixo Bento  _et al._ (2020) Bruno Valeixo Bento, Dibya Chakraborty, Susha L. Parameswaran, and Ivonne Zavala, “Dark Energy in String Theory,” PoS CORFU2019, 123 (2020), arXiv:2005.10168 [hep-th] .
  * Brinkmann  _et al._ (2022) Max Brinkmann, Michele Cicoli, Giuseppe Dibitetto, and Francisco G. Pedro, “Stringy multifield quintessence and the Swampland,” (2022), arXiv:2206.10649 [hep-th] .
  * Damour and Polyakov (1994) T. Damour and Alexander M. Polyakov, “The String dilaton and a least coupling principle,” Nucl. Phys. B 423, 532–558 (1994), arXiv:hep-th/9401069 .
  * Brax  _et al._ (2010) Philippe Brax, Carsten van de Bruck, Anne-Christine Davis, and Douglas Shaw, “The Dilaton and Modified Gravity,” Phys. Rev. D 82, 063519 (2010), arXiv:1005.3735 [astro-ph.CO] .
  * Amendola (2004) Luca Amendola, “Linear and non-linear perturbations in dark energy models,” Phys. Rev. D 69, 103524 (2004), arXiv:astro-ph/0311175 .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Lesgourgues (2011) Julien Lesgourgues, “The cosmic linear anisotropy solving system (class) i: Overview,” (2011), arXiv:1104.2932 [astro-ph.IM] .
  * Blas  _et al._ (2011) Diego Blas, Julien Lesgourgues, and Thomas Tram, “The cosmic linear anisotropy solving system (CLASS). part II: Approximation schemes,” Journal of Cosmology and Astroparticle Physics 2011, 034–034 (2011).
  * Brinckmann and Lesgourgues (2018) Thejs Brinckmann and Julien Lesgourgues, “MontePython 3: boosted MCMC sampler and other features,” (2018), arXiv:1804.07261 [astro-ph.CO] .
  * Audren  _et al._ (2013) Benjamin Audren, Julien Lesgourgues, Karim Benabed, and Simon Prunet, “Conservative Constraints on Early Cosmology: an illustration of the Monte Python cosmological parameter inference code,” JCAP 1302, 001 (2013), arXiv:1210.7183 [astro-ph.CO] .
  * Lewis (2019) Antony Lewis, “GetDist: a Python package for analysing Monte Carlo samples,”  (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Alam  _et al._ (2017) S. Alam _et al._ , “The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample,” Monthly Notices of the Royal Astronomical Society 470, 2617–2652 (2017), https://academic.oup.com/mnras/article-pdf/470/3/2617/18315003/stx721.pdf .
  * Beutler  _et al._ (2011) Florian Beutler, Chris Blake, Matthew Colless, D. Heath Jones, Lister Staveley-Smith, Lachlan Campbell, Quentin Parker, Will Saunders, and Fred Watson, “The 6dF Galaxy Survey: baryon acoustic oscillations and the local Hubble constant,” Monthly Notices of the Royal Astronomical Society 416, 3017–3032 (2011), https://academic.oup.com/mnras/article-pdf/416/4/3017/2985042/mnras0416-3017.pdf .
  * Ross  _et al._ (2015) Ashley J. Ross, Lado Samushia, Cullan Howlett, Will J. Percival, Angela Burden, and Marc Manera, “The clustering of the SDSS DR7 main Galaxy sample – I. A 4 per cent distance measure at z = 0.15,” Monthly Notices of the Royal Astronomical Society 449, 835–847 (2015), https://academic.oup.com/mnras/article-pdf/449/1/835/13767551/stv154.pdf .
  * Scolnic  _et al._ (2018) D. M. Scolnic _et al._ (Pan-STARRS1), “The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample,” Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Sagredo  _et al._ (2018) Bryan Sagredo, Savvas Nesseris, and Domenico Sapone, “Internal Robustness of Growth Rate data,” Phys. Rev. D 98, 083543 (2018), arXiv:1806.10822 [astro-ph.CO] .
  * Arjona  _et al._ (2020) Rubén Arjona, Juan García-Bellido, and Savvas Nesseris, “Cosmological constraints on non-adiabatic dark energy perturbations,” (2020), arXiv:2006.01762 [astro-ph.CO] .
  * Moresco  _et al._ (2016) Michele Moresco, Lucia Pozzetti, Andrea Cimatti, Raul Jimenez, Claudia Maraston, Licia Verde, Daniel Thomas, Annalisa Citro, Rita Tojeiro, and David Wilkinson, “A 6% measurement of the hubble parameter at $z\sim 0.45$: direct evidence of the epoch of cosmic re-acceleration,” Journal of Cosmology and Astroparticle Physics 2016, 014–014 (2016).
  * Gelman and Rubin (1992) Andrew Gelman and Donald B. Rubin, “Inference from Iterative Simulation Using Multiple Sequences,” Statistical Science 7, 457 – 472 (1992).
  * van de Bruck and Mifsud (2018) Carsten van de Bruck and Jurgen Mifsud, “Searching for dark matter - dark energy interactions: going beyond the conformal case,” Phys. Rev. D 97, 023506 (2018), arXiv:1709.04882 [astro-ph.CO] .
  * Kesden and Kamionkowski (2006) Michael Kesden and Marc Kamionkowski, “Tidal Tails Test the Equivalence Principle in the Dark Sector,” Phys. Rev. D 74, 083007 (2006), arXiv:astro-ph/0608095 .


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

Joshua M. Lewis1,2 and Lincoln D. Carr1,2,3 1Quantum Engineering Program,
Colorado School of Mines, Golden, CO 80401, U.S.A.  2Department of Physics,
Colorado School of Mines, Golden, CO 80401, U.S.A.  3Department of Applied
Mathematics and Statistics, Colorado School of Mines, Golden, CO 80401, U.S.A.
<EMAIL_ADDRESS>
###### Abstract

Fractional evolution equations lack generally accessible and well-converged
codes excepting anomalous diffusion. A particular equation of strong interest
to the growing intersection of applied mathematics and quantum information
science and technology is the fractional Schrödinger equation, which describes
sub-and super-dispersive behavior of quantum wavefunctions induced by
multiscale media. We derive a computationally efficient sixth-order split-step
numerical method to converge the eigenfunctions of the FSE to arbitrary
numerical precision for arbitrary fractional order derivative. We demonstrate
applications of this code to machine precision for classic quantum problems
such as the finite well and harmonic oscillator, which take surprising twists
due to the non-local nature of the fractional derivative. For example, the
evanescent wave tails in the finite well take a Mittag-Leffer-like form which
decay much slower than the well-known exponential from integer-order
derivative wave theories, enhancing penetration into the barrier and therefore
quantum tunneling rates. We call this effect _fractionally enhanced quantum
tunneling_. This work includes an open source code for communities from
quantum experimentalists to applied mathematicians to easily and efficiently
explore the solutions of the fractional Schrödinger equation in a wide variety
of practical potentials for potential realization in quantum tunneling
enhancement and other quantum applications.

##  1 Fractional differential equations in context

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.

##  2 Overview of the Fractional Schrödinger Equation

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.

##  3 The Mittag-Leffler Function

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.

(a)

(b)

Figure 1: _Comparison of the Mittag-Leffler function to the Gaussian._ (a):
Core behavior of the Mittag-Leffler function. (b): Logarithmic scale
emphasizing the tails.

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.

##  4 Fourier Split Step: The Fractional Adaptation

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.

Step | $a_{n}$ | $b_{n}$
---|---|---
1 |  | 0.0584500187773306
---
+0.0217141273080301i
| 0.116900037554661
---
+0.0434282546160603i
2 |  | 0.123229569418374
---
-0.0402806787860161i
| 0.129559101282088
---
-0.123989612188092i
3 |  | 0.158045797047111
---
-0.0604410907390099i
| 0.186532492812133
---
+0.00310743071007267i
4 |  | 0.160274614757183
---
+0.0790076422169959i
| 0.134016736702233
---
+0.154907853723919i
5 | $a_{4}$ | $b_{3}$
6 | $a_{3}$ | $b_{2}$
7 | $a_{2}$ | $b_{1}$
8 | $a_{1}$ | 0
Table 1: The splitting coefficients used in the FSE algorithm to obtain a
sixth order approximation.

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].

##  5 Demonstration of FSE Code Precision: Particle on a Ring

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.

Figure 2: _Point-wise error associated with the 1st excited state of the
particle on a ring system across fractional orders._

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.

Figure 3: _Maximum point-wise error across many eigenstates of the particle on
a ring system for a few fractional orders._

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

Figure 4: _Error in generated energies when compared to analytical energy._

##  6 Fractional Quantum Mechanics Eigensystems

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.

###  6.1 Quantum Harmonic Oscillator

The quantum harmonic oscillator problem is a somewhat famous problem within
the study of quantum mechanics as the potential in the system well
approximates many physically important systems and eigenstates produced lead
to the construction of the coherent state and describe many important physical
phenomena. The system be defined with the following Hamiltonian.

$\hat{H}=-\frac{1}{2}\partial_{x}^{\alpha}+\frac{1}{2}x^{2}$ (20)

Contrary to the prior system the eigenfunctions are no longer invariant of
fractional order. We compare the ground state in Fig. 5(a) across different
orders to see this. While these different ground states appear to be of the
same structure with slight variations, the greatest difference is hidden from
sight in this plot. To uncover the main difference between these ground states
we will analyze the log plot in Fig. 5(b).

The tails of the ground state are greatly altered by changing the fractional
order. We see that the exponential decay from the integer order case has been
replaced by Mittag-Leffler like decay in the tails. The node in the order 2.2
case is a clear hint of the Mittag-Leffler like nature of this eigenstate. One
may refer back to Fig. 1 to see the very similar behavior of the tails in the
classically forbidden region.

(a) __

(b) __

(c) __

Figure 5: _Eigenfunctions of the quantum harmonic oscillator potential across
fractional orders: (a) Ground states depicted across fractional orders less
than, equal to, and greater than 2. (b) Logarithmic view emphasizing tail
visibility. (c) Close-up at the transition into the classically forbidden
region, highlighting the distinct Mittag-Leffler-like tails and node in the
higher order eigenstate._

We may then recover the energies associated with each eigenfunction by again
evolving the state for a short time in imaginary time and extracting the
energy from the decay rate.

Eigenfunction | $\alpha=1.8$ | $\alpha=2.0$ | $\alpha=2.2$
---|---|---|---
0 | 0.4994984133 | 0.4999999999 | 0.5012687387
1 | 1.4452024402 | 1.5000000003 | 1.5508889832
2 | 2.3415247398 | 2.4999999995 | 2.6507789545
3 | 3.2229118004 | 3.4999999991 | 3.7689276297
4 | 4.0884863294 | 4.4999999983 | 4.9044885301
Table 2: _Energies for the first five eigenfunctions associated with the
quantum harmonic oscillator of orders 1.8, 2.0, and 2.2_

One may notice that the linear eigenspectrum afforded to the system with
integer order is entirely lost when the fractional order is changed. Thus,
changing the fractional order of the quantum harmonic oscillator results in an
eigenspectrum that does not have a revival time nor does it have a similar
coherent state. We may then plot the spectrum of eigenvalues across many
orders to see a clearer trend of how the eigenspectrum changes upon a changing
fractional order.

Figure 6: _Spectrum of energies associated with eigenfunctions of the
fractional quantum harmonic oscillator across many orders_

###  6.2 Finite Well

To examine the tails in the classically forbidden region and the compaction or
expansion associated with the spectrum of eigenenergies upon a changing
fractional order, we now investigate the finite well problem which may be
defined with the following Hamiltonian.

$\hat{H}=-\frac{1}{2}\partial_{x}^{\alpha}+V(x)$ (21)

$V(x)=\begin{cases}V_{0}&\text{otherwise}\\\ 0&\text{if }|x|<1\end{cases}$
(22)

We will choose to set $V_{0}$ to 100 in this analysis such that a sufficient
number of bound eigenstates may be generated within the well. We observe the
ground state across three different fractional orders such that adequate
comparisons may be made.

Observation of Fig. 7(a) shows that the ground state across fractional orders
appears to be changed in only minute and insignificant ways. However, once
more, the true difference lies within the tails within the classically
forbidden region as is seen in Fig. 7(b).

(a) __

(b) __

(c) __

Figure 7: _Same analysis as in Fig. 5._

We may immediately note again that the tails of these eigenstates are no
longer well characterized by exponential decay. They are in fact Mittag-
Leffler like decay patterns similar to what was obtained in the fractional
quantum harmonic oscillator and what was shown in Fig. 1. Given this drastic
change to the tails within the classically forbidden region, we may
hypothesize that there is an enhancement to quantum tunneling due to the
longer range nature of confined states.

Furthermore, we may examine the number of bound states within this system. As
is well known, eigenstates with energies less than the height of the finite
well are bound states and as such they are effectively quantized. We may then
plot how many eigenstates are bound across many orders.

Figure 8: _Number of bound eigenstates in the finite well system across many
fractional orders_

It is clear that as the order of the fractional derivative increases so does
the energy spacing between states; this results in a smaller number of states
that are able to exist within the finite well. Additionally as the fractional
order decreases, the number of states that may be placed in the same well
increase quite quickly.

Finally, in the analysis of finite well potentials using our Fourier split-
step method with a 6th order Suzuki-Trotter approximation, the adjustment of
the time step size plays a critical role in managing the error terms,
particularly those associated with the commutator in the Suzuki-Trotter
expansion. The non-commutative nature of the kinetic and potential energy
operators introduces significant error terms in the presence of
discontinuities, as found in steep finite wells. These error terms are
inherently linked to the commutator of these operators.

By decreasing the time step size, we directly influence the magnitude of the
commutator term in the error expression. In the Suzuki-Trotter approximation,
the error is proportional to the higher-order commutators of the Hamiltonian
components, which are functions of the time step size. Smaller time steps lead
to a reduced impact of these commutator terms, thereby diminishing the overall
error. This reduction is particularly crucial in the vicinity of the
discontinuities in the finite well, where the standard error terms would
otherwise lead to significant inaccuracies.

In the final stages of our simulation, the implementation of an even smaller
time step for a brief period serves to further alleviate the distortions
caused by these commutator terms. This meticulous calibration of the time step
ensures that the evolution of the wavefunction remains accurate, despite the
challenges posed by the non-smooth nature of the potential.

###  6.3 Double Well

As a practical example of the enhancement of quantum tunnelling we investigate
the quantum double well. A state localized in one well with less energy than
the barrier separating them are able to tunnel into the other well through the
classically forbidden zone in the middle. We may define the Hamiltonian of
this system as the following.

$\hat{H}=-\frac{1}{2}\partial_{x}^{\alpha}+V(x)$ (23)

$V(x)=-4x^{2}+\frac{1}{2}x^{4}+8$ (24)

While in the algorithm a large range is chosen as to make sure the
eigenfunction is not cut off in the classical forbidden zone we may plot a
close up view of the potential below. We also plot the two lowest energy
eigenstates which are a symmetric and anti-symmetric wavefunction for order
1.9 as an example.

Figure 9: _Potential for double well system along with the two lowest energy
eigenstates for fractional order 1.9 as an example_

Given these two states, we may put them into an equal superposition that
almost entirely isolates it to a specific well. In fact we may construct two
superpositions, one in which the wavefunction is in the left well and one that
is in the right well. We may write them as the following.

$\lvert\psi_{L}\rangle=\frac{1}{\sqrt{2}}(\lvert\psi_{+}\rangle-\lvert\psi_{-}\rangle)$
(25)

$\lvert\psi_{R}\rangle=\frac{1}{\sqrt{2}}(\lvert\psi_{+}\rangle+\lvert\psi_{-}\rangle)$
(26)

Where $\lvert\psi_{+}\rangle$ is the symmetric wavefunction and
$\lvert\psi_{-}\rangle$ is the anti-symmetric wavefunction. The reason that
tunnelling occurs is that the energy of the symmetric wavefunction is just
slightly different from the anti-symmetric wavefunction. Taking the
wavefunction in the left well, if we choose to evolve the system in time it
will take on the following form.

$e^{-i\hat{H}t}\lvert\psi_{L}\rangle=\frac{1}{\sqrt{2}}(e^{-iE_{0}t}\lvert\psi_{+}\rangle-e^{-iE_{1}t}\lvert\psi_{-}\rangle)$
(27)

Removing the unimportant total phase we obtain the following.

$e^{-i\hat{H}t}\lvert\psi_{L}\rangle=\frac{1}{\sqrt{2}}(\lvert\psi_{+}\rangle-e^{-it(E_{1}-E_{0})}\lvert\psi_{-}\rangle)$
(28)

$e^{-i\hat{H}\frac{\pi}{(E_{1}-E_{0})}}\lvert\psi_{L}\rangle=\lvert\psi_{R}\rangle$
(29)

Clearly after at a time of $t=\frac{\pi}{(E_{1}-E_{0})}$ the wavefunction
isolated to the left well becomes isolated to the right well as the relative
phase between these two states changes. After the same amount of time it will
then return to the left well again. This means that these states will
oscillate tunnelling between wells with the following frequency.

$f=\frac{|E_{1}-E_{0}|}{2\pi}$ (30)

We may now show the differences in energies for the first two eigenstates of
the double well system across many fractional orders.

Figure 10: _Energy difference between the two lowest energy states across many
fractional orders_

As can be seen, the energy difference for this system with a fractional order
less than 2.0 is larger and thus in those systems tunnelling occurs at a
greater frequency. Interestingly enough when one increases the fractional
order at first tunnelling frequency decreases, before rebounding and
increasing again similar to the case where one decreased the fractional order.
While this plot looks very linear, it is just slightly not with a small
negative second derivative for the data on the order of 10e-6.

##  7 Conclusion

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.

## References

  * Sun et al. [2018] HongGuang Sun, Yong Zhang, Dumitru Baleanu, Wen Chen, and YangQuan Chen.  A new collection of real world applications of fractional calculus in science and engineering.  _Communications in Nonlinear Science and Numerical Simulation_ , 64:213–231, 2018.
  * Kulish and Lage [2002] Vladimir V Kulish and José L Lage.  Application of fractional calculus to fluid mechanics.  _J. Fluids Eng._ , 124(3):803–806, 2002.
  * Wyss [1986] Walter Wyss.  The fractional diffusion equation.  _Journal of Mathematical Physics_ , 27(11):2782–2785, 1986.
  * Sokolov et al. [2002] Igor M Sokolov, Joseph Klafter, and Alexander Blumen.  Fractional kinetics.  _Physics Today_ , 55(11):48–54, 2002.
  * Lai et al. [2016] Jinxing Lai, Sheng Mao, Junling Qiu, Haobo Fan, Qian Zhang, Zhinan Hu, Jianxun Chen, et al.  Investigation progresses and applications of fractional derivative model in geotechnical engineering.  _Mathematical Problems in Engineering_ , 2016, 2016.
  * Chen et al. [2022] Wen Chen, HongGuang Sun, Xicheng Li, et al.  _Fractional derivative modeling in mechanics and engineering_.  Springer, 2022.
  * Obembe et al. [2017] Abiola D Obembe, Hasan Y Al-Yousef, M Enamul Hossain, and Sidqi A Abu-Khamsin.  Fractional derivatives and their applications in reservoir engineering problems: A review.  _Journal of Petroleum Science and Engineering_ , 157:312–327, 2017.
  * Zhang et al. [2021] Qiao Zhang, Yuxin Sun, and Jialing Yang.  Bio-heat transfer analysis based on fractional derivative and memory-dependent derivative heat conduction models.  _Case Studies in Thermal Engineering_ , 27:101211, 2021\.
  * Ionescu et al. [2017] C Ionescu, A Lopes, Dana Copot, JA Tenreiro Machado, and Jason HT Bates.  The role of fractional calculus in modeling biological phenomena: A review.  _Communications in Nonlinear Science and Numerical Simulation_ , 51:141–159, 2017.
  * Zhokh et al. [2017] AA Zhokh, AI Trypolskyi, and PE Strizhak.  An investigation of anomalous time-fractional diffusion of isopropyl alcohol in mesoporous silica.  _International Journal of Heat and Mass Transfer_ , 104:493–502, 2017.
  * Zhokh and Strizhak [2017] Alexey A Zhokh and Peter E Strizhak.  Experimental verification of the time-fractional diffusion of methanol in silica.  _J. Appl. Nonlinear Dyn_ , 6(2):135, 2017.
  * Chen and Raghavan [2015] C Chen and R Raghavan.  Transient flow in a linear reservoir for space–time fractional diffusion.  _Journal of Petroleum science and engineering_ , 128:194–202, 2015.
  * Raghavan [2011] R Raghavan.  Fractional derivatives: Application to transient flow.  _Journal of Petroleum Science and Engineering_ , 80(1):7–13, 2011.
  * Mainardi et al. [2007] Francesco Mainardi, Gianni Pagnini, and Rudolf Gorenflo.  Mellin transform and subordination laws in fractional diffusion processes.  _arXiv preprint math/0702133_ , 2007.
  * Abdel-Rady et al. [2015] Ahmed Safwat Abdel-Rady, Saad Zagloul Rida, Anas Ahmed Mohamed Arafa, Hamdy Ragab Abedl-Rahim, et al.  Natural transform for solving fractional models.  _Journal of Applied Mathematics and Physics_ , 3(12):1633, 2015.
  * Leforestier et al. [1991] Claude Leforestier, RH Bisseling, Charly Cerjan, MD Feit, Rich Friesner, A Guldberg, A Hammerich, G Jolicard, W Karrlein, H-D Meyer, et al.  A comparison of different propagation schemes for the time dependent schrödinger equation.  _Journal of Computational Physics_ , 94(1):59–80, 1991.
  * Klauder [2003] John R Klauder.  The feynman path integral: An historical slice.  In _A Garden of Quanta: Essays in Honor of Hiroshi Ezawa_ , pages 55–76. World Scientific, 2003.
  * Laskin [2000a] Nick Laskin.  Fractional quantum mechanics.  _Physical Review E_ , 62(3):3135, 2000a.
  * Sun et al. [2010] HongGuang Sun, Wen Chen, Changpin Li, and YangQuan Chen.  Fractional differential models for anomalous diffusion.  _Physica A: statistical mechanics and its applications_ , 389(14):2719–2724, 2010.
  * Laskin [2000b] Nick Laskin.  Fractals and quantum mechanics.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 10(4):780–790, 2000b.
  * Laskin [2002] Nick Laskin.  Fractional schrödinger equation.  _Physical Review E_ , 66(5):056108, 2002.
  * Laskin [2010] Nick Laskin.  Principles of fractional quantum mechanics.  _arXiv preprint arXiv:1009.5533_ , 2010.
  * Jeng et al. [2010] M Jeng, S-L-Y Xu, Eli Hawkins, and JM Schwarz.  On the nonlocality of the fractional schrödinger equation.  _Journal of mathematical physics_ , 51(6), 2010.
  * Akhiezer and Glazman [2013] Naum Il’ich Akhiezer and Izrail Markovich Glazman.  _Theory of linear operators in Hilbert space_.  Courier Corporation, 2013.
  * Kilbas et al. [2004] Anatoly A Kilbas, Megumi Saigo, and Ram K Saxena.  Generalized mittag-leffler function and generalized fractional calculus operators.  _Integral Transforms and Special Functions_ , 15(1):31–49, 2004.
  * Rahman et al. [2017] Gauhar Rahman, Dumitru Baleanu, Maysaa Al Qurashi, Sunil Dutt Purohit, Shahid Mubeen, and Muhammad Arshad.  The extended mittag-leffler function via fractional calculus.  _J. Nonlinear Sci. Appl_ , 10(8):4244–4253, 2017\.
  * Peng and Li [2010] Jigen Peng and Kexue Li.  A note on property of the mittag-leffler function.  _Journal of Mathematical Analysis and Applications_ , 370(2):635–638, 2010.
  * Blanes and Casas [2005] Sergio Blanes and Fernando Casas.  On the necessity of negative coefficients for operator splitting schemes of order higher than two.  _Applied Numerical Mathematics_ , 54(1):23–37, 2005.
  * Lewis [2024] Joshua Lewis.  Fractionaleigen: Source code for fourier split-step simulations.  https://github.com/JoshLewisGH/FractionalEigen, 2024.