Fractional Path Integral Monte Carlo

# Fractional Path Integral Monte Carlo

## Abstract

Fractional derivatives are nonlocal differential operators of real order that often appear in models of anomalous diffusion and a variety of nonlocal phenomena. Recently, a version of the Schrödinger Equation containing a fractional Laplacian has been proposed. In this work, we develop a Fractional Path Integral Monte Carlo algorithm that can be used to study the finite temperature behavior of the time-independent Fractional Schrödinger Equation for a variety of potentials. In so doing, we derive an analytic form for the finite temperature fractional free particle density matrix and demonstrate how it can be sampled to acquire new sets of particle positions. We employ this algorithm to simulate both the free particle and He (Aziz) Hamiltonians. We find that the fractional Laplacian strongly encourages particle delocalization, even in the presence of interactions, suggesting that fractional Hamiltonians may manifest atypical forms of condensation. Our work opens the door to studying fractional Hamiltonians with arbitrarily complex potentials that escape analytical solutions.

Path Integrals, Fractional Schrödinger Equation, Fractional Laplacian, Anomalous Diffusion, Monte Carlo

## I Introduction

Fractional calculus is the study and use of differential operators of real order. Such operators are nonlocal except for the case of nonnegative integer order. Space- and time-fractional differential equations have emerged as the correct way to model anomalous diffusion, Klafter and Sokolov (2005); Meerschaert and Sikorskii (2011) nonlocal phenomena,Lazopoulos (2006) memory effects,Du, Wang, and Hu (2013) and many non-equilibrium Vlahos et al. (2008) and turbulentMajda and Kramer (1999); Richardson (1926) systems.

Anomalous diffusion is any diffusion in which the mean square displacement is not linear in time, i.e.,

 ⟨x2⟩∼t2α, (1)

where , in contrast to the case, which corresponds to normal/Brownian diffusion. If , the diffusion is superdiffusive, while if , the diffusion is subdiffusive. One can think of superdiffusion for as a regime between normal diffusion and ballistic motion (for which , implying that the average displacement is proportional to time). A simple example of superdiffusion is a plasma whose particles would ordinarily be expected to undergo diffusion because of their fluid nature, but instead exhibit superdiffusive behavior because of the presence of a magnetic field that favors ballistic motion.Corkum (1973); Vlahos et al. (2008) On the other hand, traps (such as in porous mediaBenson et al. (2001)) or formations of eddies in a fluidVlahos et al. (2008); Solomon, Weeks, and Swinney (1994); Weeks, Urbach, and Swinney (1996) may introduce waiting times that slow down the fluid’s particles, pushing their motion toward subdiffusion. In general, both tendencies can occur, resulting in a unique process with “competition” between super- and sub-diffusive overall behavior.

Anomalous diffusion differs significantly from normal diffusion, in more ways than can be expected from the scaling relation given by Relation (1). Anomalous diffusion is often characterized by heavy tails, infinite variance, and inhomogeneity in space and time. These are characteristics of Continuous Time Random Walks (CTRWs), a broad family of anomalous diffusions which model both waitings times and long jumps. Meerschaert and Scheffler (2004); Meerschaert and Sikorskii (2011) The simplest subfamily of CTRWs is isotropic -stable Levy motion.

Just as the standard Laplacian in is the generator of Brownian motion, a random walk in which the jumps can be drawn from a normal distribution, for , the fractional Laplacian in is the generator of isotropic -stable Levy motion.Meerschaert and Sikorskii (2011) In this motion, jumps are drawn from an -stable distribution, which exhibits power-law tails. The isotropic -stable distribution is defined by a position parameter and a scale parameter, which for the normal case, reduces to the mean and variance, respectively.

The field of physics where anomalous diffusion is most widespread is perhaps plasma physics, where anomalous diffusion has been observed for decades.Yoshikawa and Rose (1962); Geissler (1968) This has prompted theoretical discussions on how the equations of magnetohydrodynamics can drive anomalous diffusionDrummond and Rosenbluth (1962); Corkum (1973); Zaslavsky (2002) and applications of CTRWs to plasma physics.Balescu (1995) Anomalous diffusion has moreover been observedKapitulnik and Deutscher (1983) in percolation clusters, which has since been explained theoretically. Gefen, Aharony, and Alexander (1983) Inspired by this finding, the role of anomalous diffusion in disordered superconductors (and other disordered media) has also been elucidated.Bouchaud and Georges (1990)

Recently, a nonlinear fractional Schrödinger equation has been rigorously derived from a mathematical physics perspective as the continuum limit of a system of discrete nonlinear Schrödinger equations with increasingly long-range interactions.Kirkpatrick, Lenzmann, and Staffilani (2013) LaskinLaskin (2000, 2002, 2008, 2010) first developed fractional quantum mechanics, following the early considerations of Kac.Kac (1951) Starting from the Feynman path-integral formulation of quantum mechanics, Laskin replaced integrals over Brownian paths with integrals over -stable Levy paths. He showed that this amounts to replacing the standard Laplacian by the (Riesz)1 fractional Laplacian in the Schrödinger Equation, leading to

 iℏ∂ψ(→r,t)∂t=−ℏDαΔα/2ψ(→r,t)+V(→r)ψ(→r,t). (2)

The equation is called the space-fractional Schrödinger equation. Laskin showed that this equation enjoys the same basic properties that allow the usual Schrödinger equation to serve as a foundation for quantum mechanics; for example, since is self-adjoint, observables are real valued.Laskin (2010) To date, solutions to the fractional Schrödinger Equation have been derived for a variety of analytically tractable Hamiltonians, including the free particle, square well, and particle in a box Hamiltonians.Guo and Xu (2006); Dong and Xu (2008a); Secchi (2013); Felmer, Quaas, and Tan (2012) Solutions to Hamiltonians containing more realistic potentials have so far been few and far between, limiting the study of potentially rich fractional mathematics and physics.

In this paper, we introduce a new fractional Path Integral Monte Carlo method capable of modeling the finite temperature physics of a wide variety of space-fractional Hamiltonians. The method computes finite temperature observables using the Path Integral Monte Carlo (PIMC) algorithm popularly used to characterize boson and Boltzmannon phase diagrams Ceperley (1995); Ceperley and Pollock (1986); Pollock and Ceperley (1984) to sample the many body fractional partition function. Key to being able to use the PIMC algorithm in this context is being able to express and sample the free finite temperature fractional density matrix. We therefore derive an analytical form for this fractional density matrix and show how it can be both computed and sampled.

Interestingly, we find that the probability distribution function that represents the fractional density matrix manifests dramatically different asymptotic behaviors depending upon the fractional exponent present in the Hamiltonian. Whereas for fractional exponent the distribution reduces to the usual Gaussian case, for , the distribution possesses heavy tails, while for , the distribution becomes negative. Setting aside the case, which is typically excluded2 in discussions of the fractional laplacian, in the remainder of this work, we illustrate how our algorithm may be used to compute finite temperature observables including energies and radii of gyration for two illustrative Hamiltonians containing fractional Laplacians: 1) the free particle Hamiltonian; and 2) the He Hamiltonian using the AzizAziz et al. (1979) potential. We end with a discussion of the general applicability of our method and the intriguing fractional physics it may reveal.

## Ii Theory and Algorithms

### ii.1 Overview of the Path Integral Monte Carlo Method

Since our method may be viewed as a generalization of the conventional Path Integral Monte Carlo algorithm, we begin this section with a review of PIMC. As described in significantly more detail by Ceperley,Ceperley (1995) the finite temperature equilibrium expectation values of a system observable, , in the position basis may be expressed as

 ⟨^O⟩=Z−1∫d→Rd→R′ρ(→R,→R′;β)⟨→R|^O|→R′⟩, (3)

where denotes the system’s partition function

 Z=∫d→Rρ(→R,→R;β) (4)

at inverse temperature . In the above, is the position vector of all particles’ coordinates, where represents the position of the th particle. One of the key quantities present in both Equations (3) and (4) is , the finite temperature density matrix. In position space, the density matrix may be expanded into

 ρ(→R,→R′;β) = ⟨→R|e−β^H|→R′⟩ (5) = ∑iϕ∗i(→R)ϕi(→R′)e−βEi,

where is the Hamiltonian that describes the system, and and respectively denote its eigenvalues (energies) and eigenfunctions (wave functions). Although obtaining an exact expression for the density matrix at an arbitrary temperature would require diagonalizing the many body Hamiltonian, a task generally beyond our current computational capabilities, the density matrix may be simplified into more tractable expressions via Trotter factorization. In general, the convolution of two density matrices remains a density matrix such that

 ρ(→R1,→R3;β1+β2)=∫d→R2ρ(→R1,→R2;β1)ρ(→R2,→R3;β2), (6)

which implies that the expression for the density matrix at inverse temperature, , may be expressed as the convolution of density matrices at inverse temperature

 ρ(→R0,→RM;β) = ∫...∫d→R1d→R2...d→RM−1ρ(→R0,→R1;τ) (7) ρ(→R1,→R2;τ)...ρ(→RM−1,→RM;τ).

The different -particle position vectors represent the particle positions at different so-called “imaginary times.” Because is the same as in the expression for the partition function, according to the quantum-classical isomorphism, the particles may be thought of as polymers constructed of different beads linked together. If is sufficiently small, assuming the system Hamiltonian may be split into kinetic () and potential pieces (), , the Suzuki-Trotter factorizationTrotter (1959); Kato (1974) may be used to break up the exponential of the Hamiltonian into corresponding kinetic and potential exponentials

 e−τ^H=e−τ(^K+^V)≈e−τ^Ke−τ^V (8)

such that the position-space density matrices may also be broken into

 ρ(→R0,→R2;τ)≡∫d→R1⟨→R0|e−τ^K|→R1⟩⟨→R1|e−τ^V|→R2⟩. (9)

Because potential operators are generally diagonal in the position representation, evaluating the position density matrix is trivial

 ⟨→R1|e−τ^V|→R2⟩=e−τ^V(→R1)δ(→R2−→R1). (10)

The evaluation of the kinetic density matrix is less trivial and requires finding the eigenvalues and eigenfunctions of . Assuming periodic boundary conditions and the , non-fractional form for the kinetic operator

 ^K=−ℏ22m∇2=−λ∇2, (11)

in which denotes the Laplacian, denotes the differential operator, and , the eigenvalues of the three-dimensional kinetic operator are and the eigenfunctions are , with . The kinetic density matrix may therefore be expressed as

 ⟨→R0|e−τ^K|→R1⟩ = ∑→nL−3Ne−τλK2n−i→Kn(→R0−→R1) (12) = (4πλτ)−3N/2e−(→R0−→R1)24λτ.

As discussed in more detail in Section II.6.1, the finite temperature partition function may then be sampled via Monte Carlo by sampling the Gaussian kinetic density matrix probability distribution function for new particle coordinates and then using the potential density matrices to accept/reject those coordinates according to their potential energy contributions. As a prelude to the subsequent discussion, it should be noted that Equation (12), and the eigenvalues and eigenfunctions used to evaluate it, is only valid for the Laplacian. If the differential operator is of a different power, the Suzuki-Trotter factorization is still valid (see Appendix II), but new expressions for the density matrix are required. These expressions are derived and discussed below.

### ii.2 The Fractional Schrödinger Equation

As first presented by Laskin,Laskin (2000, 2002, 2008, 2010) the fractional Schrödinger Equation is a generalization of the Schrödinger Equation in which the Hamiltonian is generalized to a fractional form, , that contains a fractional Laplacian, . In one dimension, , may be expressed as

 ^Hα=−DαℏαΔα/2+^V(x). (13)

In the above, denotes the potential and denotes the constant that falls in front of the Laplacian. In general, must have SI units of

 [Dα]=m2−αkg1−αs2−α. (14)

The fractional exponent, , may in principle assume any positive value, including non-integer, fractional values, as its name implies.Laskin (2000, 2002, 2008, 2010) For reference, the fractional Hamiltonian reduces to the conventional case when . Please also note that in the next few sections we consider the one-dimensional version of the Hamiltonian to simplify our discussion of the derivation of the fractional density matrix; we will return to a discussion of the full three-dimensional form in later sections.

Most previous work on the fractional Schrödinger Equation has explored the properties of fractional Hamiltonians with potentials that are analytically soluble, including the delta potential,de Oliveira, Costa, and Jr. (2010), the fractional hydrogen atom,Laskin (2002) the fractional oscillator,Laskin (2002) and the fractional square well.Guo and Xu (2006); Dong and Xu (2008a); Secchi (2013); Felmer, Quaas, and Tan (2012) In addition to the fact that these models manifest atypical physics, such as atypical energy spectra, they are also of intrinsic mathematical interest. Considering that that the ordinary Schrödinger equation is already extremely challenging to solve, the inclusion of a (nonlocal) fractional derivative makes equation (2) even less tractable. More recently, it has been suggested that fractional Hamiltonians may naturally arise from strong many body interactions, which prompted a mean field exploration of fractional Bose-Einstein condensation.Kleinert (2012) Nevertheless, no robust algorithm for examining the properties of the fractional Schrödinger Equation for arbitrary potentials has yet been proposed.

### ii.3 Spectrum of the Fractional Laplacian

In order to generalize the Path Integral Monte Carlo algorithm discussed in Section II.1 to fractional Hamiltonians, expressions for the fractional density matrix must be obtained. As discussed above, this requires finding an analytical form for the kinetic (non-interacting) density matrix that is also amenable to sampling. This can be achieved by finding the eigenvalues and eigenfunctions of the fractional Laplacian. The eigenvalue spectrum of the Riesz fractional Laplacian in one dimension is given by

 Δα/2cos(Cx) = −|C|αcos(Cx) (15) Δα/2sin(Cx) = −|C|αsin(Cx) (16)

for any in the unbounded case.Hermann (2014) In fact, this can be viewed as the definition of the real power of the Laplacian operator via the spectral theorem.Rudin (1991)

The standard boundary conditions used in condensed phase simulations are periodic boundary conditions. The introduction of periodic boundary conditions restricts the continuous family of eigenfunctions to a discrete family of eigenfunctions that satisfy the required periodicity.3 For a periodic box of length , , where , the related eigenvalues of the Laplacian are therefore

 λ=−|C|α=−(2πkL)α (17)

and the related eigenfunctions are and . As in the standard derivation for the particle in a box Hamiltonian,Krauth (2006) the sums and differences of these trigonometric functions may be combined and normalized to yield eigenfunctions of the form

 ϕper,Lk(x)=√1Le−iCx. (18)

Inserting these expressions into Equation (13), the final Hamiltonian eigenvalues (energies) may be obtained

 −DαℏαΔα/2ϕper,Lk(x) = Dαℏα|C|α√1Le−iCx (19) = Eper,Lk√1Le−iCx.

Thus,

 Eper,Lk=Dαℏα|C|α=Dα∣∣∣2πkℏL∣∣∣α. (20)

Consequently, the energies in the fractional periodic setting are similar in form to the , non-fractional case, except that the frequencies are now raised to a fractional power instead of being squared.

### ii.4 Derivation of the Fractional Density Matrix

Following Ceperley and Krauth,Ceperley (1995); Krauth (2006) the kinetic density matrix in one dimension may be obtained along the same lines as in Equation (12) by summing over the Hamiltonian’s eigenvalues and eigenfunctions

 ρper,L(x,x′,β) = ∞∑k=−∞ϕper,Lk(x)e−βEper,Lkϕ∗,per,Lk(x′) (21) = 1L∞∑k=−∞eiC(x′−x)e−βDαℏα|C|α

As usual, these sums may be turned into integrals by substituting the implied in the summation with and letting the box size, , go to infinity

 ρper,L(x,x′,β)=limL→∞1L∑C=...,−2πL,0,2πL,...(ΔCL2π)eiC(x−x′)e−βDαℏα|C|α=12π∫∞−∞dCeiC(x−x′)e−βDαℏα|C|α. (22)

If , this integral boils down to the usual Gaussian case in which completing the square yields the final line of Equation (12).

For more general values of , it is impossible to compute this integral in closed form, but it can be represented as a special function defined by a power series. The integral may be identified as the characteristic function of the generalized normal random variable . This is the random variable given by the probability density function

 fκ(μ,σ;x)=ξ(κ)=κ2σΓ(1/κ)e−∣∣x−μσ∣∣κ. (23)

The characteristic function of for was first computed explicitly in terms of special functions by Pogany and Nadarajah in 2006.PogÃÂ¡ny and Nadarajah (2010) It may be expressed as

 E{eitXGN}=√πeitμΓ(1/κ)Ψ11[(1/κ,2/κ);(1/2,1);−(σt)24], (24)

where denotes the Fox-Wright generalized hypergeometric function. This Fox-Wright function is a generalization of the generalized hypergeometric function, and is a special case of the Fox -function. It is specified by an arbitrary number, , of parameters, and is defined as the series

 Ψpq[(α1,A1),...,(αp,Ap);(β1,B1),...,(βq,Bq);z]=∞∑n=0∏pj=1Γ(αj+Ajn)∏qj=1Γ(βj+Bjn)znn!. (25)

Comparing the generalized normal distribution in the expression for the density matrix given by Equation (22) with the definition of the generalized normal distribution given by Equation (23), we have , , and . Making these substitutions, we obtain

 ρper,L (x,x′,β) (26) = 12π∫∞−∞dCeiC(x−x′)e−βDαℏα|C|α = 12π2Γ(1/α)β1/αD1/ααℏα√πei(x−x′)μΓ(1/α) ×Ψ11[(1/α,2/α);(1/2,1);−(x−x′)24β2/αD2/ααℏ2] = 1√πβ1/αD1/ασℏα ×∞∑n=0Γ(1/α+(2/α)n)Γ(1/2+n)[−(x−x′)24β2/αD2/ααℏ2]nn!.

The generalized expression for the fractional kinetic density matrix is therefore proportional to the Fox-Wright function of argument . Similar representations were obtained by LaskinLaskin (2010) for the free-particle density matrix on without boundary conditions. Laskin used the more general Fox -function in his representation, which for the given indices reduces to the Fox-Wright function shown here. In another context, similar results were derived by Mainardi and Pagnini,Mainardi and Pagnini (2007) where the integral form of (26) arises as the fundamental solution to the space-fractional diffusion equation.

At this point, an efficient way to sample this Fox-Wright density is required to properly sample the kinetic density matrix.

### ii.5 Computing and Sampling the Fractional Density Matrix

The Fox-Wright density given by Equation (26) is interpreted as a probability density function (PDF) of the jump distance, , for the random walkers. For , the Fox-Wright density reduces to a Gaussian, which is typically sampled using the Box-Muller transformation.Scott (2011)

For simplicity and portability, we have used discrete inverse transform sampling to sample the distribution (26). The method is entirely sufficient for the initial computations presented here. It is both fast and the most versatile method; since it does not rely on any special transformations, it can be used to sample any random variable for which the probability density function can be computed reliably. Our method is summarized in Figure 2 and described in full detail in this section.

The parameters , , and are not essential to the sampling methodology. We begin by reducing them to unity. If the density given in Equation (22) is denoted by , then

 ρ1,1,1(x−x′)=12π∫∞−∞eiC(x−x′)e−|C|αdC. (27)

By rescaling, may be reduced to . Letting and , so that , we get

 ρβ,Dα,ℏ(x−x′) = 12π∫∞−∞(βDαℏα)−1/αe−C′(βDαℏα)−1/α(x−x′)e−|C′|αdC′ = 12π∫∞−∞(βDαℏα)−1/αe−C′(x−x′(βDαℏα)1/α)e−|C′|αdC′ = (βDαℏα)−1/αρ1,1,1(x−x′(βDαℏα)1/α)

Therefore, in what follows, we can assume without loss of generality that .

As mentioned, the PDF can be sampled using inverse transform sampling, which requires that the related cumulative density function (CDF) be inverted. The general result is that, if a random variable has CDF , then in distribution,Rubinstein (1981) where denotes the uniform random variable on . In our case, the CDF of the random variable given by does not admit analytic inversion, so various numerical methods can be used. Options include binary search on a finely-discretized function table or a Newton method. Since PIMC simulations require samples of each of the millions of times a bead is displaced, the numerical inversion must be carried out as efficiently as possible. Any method for inverting the CDF will require accurate computation of the CDF, and therefore the PDF .

There are two direct ways to compute the PDF of : direct quadrature of the integral in the first line of Equation (26) or truncation of the Taylor series in the last line of Equation (26). The truncated Taylor series is a convenient method when the argument of the distribution is small. For larger arguments, the series (which is an expansion at ) requires increasingly higher-degree terms to be accurate, but this leads quickly to numerical blow-up. The combination of high-degree monomials, large , and gamma functions and factorials with large arguments results in overflow and precision issues. One work around is the following numerical “order of operations” trick:

 xnn!=(x1)(x2)...(xn−1)(xn). (29)

This prevents the overflow arising from computing the very large numbers and separately before dividing by pairing each large factor with a large factor , . With this trick, one can use the Taylor series to compute the Fox-Wright function accurately and quickly for arguments up to about 10, using 100 terms in the series. Beyond that point, numerical issues become more difficult to resolve and the computational cost too high.

Thus, for larger arguments, direct quadrature becomes more favorable. Due to the rapid decay of , the truncated integral

 12π∫L−LeiC(x−x′)e−|C|αdC=12π∫L−Lcos[C(x−x′)]e−|C|αdC (30)

converges rapidly in . At , the integrand is far below machine precision. Note that in Equation (30), the complex exponential has been replaced by a cosine because the related sine term is odd and therefore does not contribute to the final integral. The truncated integral can then be approximated by any quadrature rule, e.g., the trapezoid rule. For large , this is appealing as there are no overflow/precision problems. The only hazard is that the frequency of the integrand is proportional to , so the quadrature must be performed with care to ensure that the result is converged in the number of quadrature points for large arguments. Quadrature is relatively costly in terms of computational time.

Our final algorithm for sampling is as follows. Treating the PDF as discrete, we compute the CDF on in intervals of 0.01. This is done using the power series for and quadrature for . The choice of at which value to transition from using the power series to using quadrature is immaterial as long as both techniques are accurate in the given region. For , where the Fox-Wright density approaches machine precision, we truncate. The resulting CDF is normalized to correct the minute mass lost by truncation and stored in memory. To generate a random sample , we use binary search on the stored table to invert the equation , where is drawn from a uniform distribution on . This method incurs high overhead when initializing the discrete Fox-Wright CDF on a fine grid, but avoids computation of the Fox-Wright density after initialization. Because the Fox-Wright CDF is computed before Monte Carlo moves are made, this approach avoids the computational bottlenecks the computation of the CDF would otherwise present.

It is not necessary to fully discretize the density matrix; in principle, one could use a Newton solver with a coarse binary search for the initial guess to invert the CDF. However, we find the computational overhead to be similar to that needed for discretization/a binary search. In the fully discrete approach, after initialization, no expensive computations of the Fox-Wright density are ever performed, while in the Newton method, the Fox-Wright density must be computed several times at new values to generate each sample. If a specific computation proves to be extremely sensitive to the heavy tails of the Fox-Wright distribution, and it is desired to avoid truncation, the more expensive Newton method may be required.

We mention that the inversion of the Fox-Wright CDF could be simplified and accelerated if a simple asymptotic expression as a power law were known. Since this function enters into a power-law regime fairly quickly for (see Figures 3 and 4), this is a reasonable hope. Any asymptotic expansion that admits closed-form, analytic inversion would make inversion of the Fox-Wright CDF trivial for many samples. Although there exists a vast literature on asymptotic expansions for such special functions, we have not found the desired expression that will serve this purpose.

Finally, we remark that inverse-transform sampling is not the only method that applies here. The CMS (Chalmers-Mallows-Stuck) transform method Chambers, Mallows, and Stuck (1976) for sampling -stable variables does not rely on the explicit computation of the PDF given by Equation (26). The CMS method can be thought of the closest approach to the Box-Muller transform for stable procsses. However, as the developers of that method note,Chambers, Mallows, and Stuck (1976) it is not clear what method is best for very large-scale Monte Carlo simulations. This is an interesting point for future survey. Regardless, as shown below, sample generation using discrete inverse transform sampling is quite fast after the CDF has been initialized, generating samples per second on a single core. Accurate tables of the CDF can be saved, so initialization is not required for future calculations. Moreover since the discrete-inverse transform method does not rely on special transforms, it can be readily applied to an arbitrary distribution. This is the main advantage of inverse-transform sampling. For example, if it was desired to use tempered distributionsMeerschaert, Sabzikar, and Chen (2015) (see section IV) rather than pure -stable distributions, it would only be required to swap the PDF subroutines. Thus, the implementation of the inverse-transform method is worthwhile for any future work on fractional PIMC.

### ii.6 Path Integral Monte Carlo Simulations Based Upon the Fractional Density Matrix

#### Performing PIMC Moves

As described in greater detail elsewhere,Boninsegni (2005) in Path Integral Monte Carlo, a random walk is performed through the space of -particle paths at the different time slices, . A random walk may be constructed by randomly displacing the positions of the different beads within the particles. While many different move possibilities exist, one of the simplest moves is to displace a subset of the beads in the -bead path of a given particle. Let the portion of a particle ’s path to be updated be denoted by , where is the number of beads in the path to be displaced. Let denote the old set of particle bead coordinates and let denote the updated (new) set of particle bead coordinates. Then, based upon Equation (4), the probability, , that this displacement will be accepted from sampling the partition function is

 P = ∏s−1j=0ρK(→r′k+j,→r′k+j+1,Δτ)∏s−1j=0ρK(→rk+j,→rk+j+1,Δτ) (31) ×e−Δτ∑s−1j=1V(→R′k+j)e−Δτ∑s−1j=1V(→Rk+j)T(→X′→→X)T(→X→→X′).

Here, represents the transition probability, which can be selected with great freedom. The most convenient choice for the transition probability is to set it equal to the product of the kinetic density matrices, so that it can be evaluated using the expressions derived above and can cancel out other contributions to . If this choice is made,

 T(→X→→X′)=s−1∏j=0ρK(→r′k+j,→r′k+j+1,Δτ), (32)

and Equation (31) reduces to

 P=e−Δτ∑s−1j=1(V(→R′k+j)−V(→Rk+j)). (33)

These equations imply that, in a practical PIMC simulation, coordinates should first be sampled from the product of the kinetic density matrices (the length of that product depends upon the number of beads to be moved) and then the move should be accepted/rejected based upon the ratio of the potential energy in the new configuration to the potential energy in the old configuration. There are many ways to sample the transition probability and construct multiple bead moves.Sprik, Klein, and Chandler (1985); Ceperley (1995); Doll, Coalson, and Freeman (1985) One of the most common ways to propose multi-bead moves in the typical, case is to use the staging algorithm,Sprik, Klein, and Chandler (1985) which samples a hierarchy of position moves based upon the fact that a product of Gaussians is a Gaussian. Because our fractional kinetic density matrices are not in general Gaussians and products of Fox-Wright functions do not yield Fox-Wright functions, the staging algorithm cannot be straightforwardly adapted to this formalism. For the purpose of this paper (we discuss alternatives capable of more efficiently sampling multi-bead and permutation moves may in our Conclusions), we therefore use single-bead and center of mass moves, despite the ergodicity problems single-bead moves may cause at low temperatures, particularly for strong interactions. For center of mass moves, all beads in each particle are moved at once, as usual. For single-bead moves, a bead on a particle is randomly selected. New unscaled ’, ’, and ’ coordinates are then independently sampled from their Fox-Wright Distributions as described in Section II.4 and rescaled by

 σα=ℏτ1/α21/αm1/α (34)

to yield coordinates of the correct dimensions. As further detailed in Appendix I, the , , and coordinates may be independently sampled because the many particle density matrices corresponding to the free fractional Hamiltonian may be re-expressed as the product of single particle density matrices in each dimension. The resulting potential energy from this move is lastly determined and used to accept/reject the proposed move.

#### Computing Observables

The fundamental purpose of performing PIMC simulations in most cases is to obtain measures of different systems’ finite temperature observables, such as their energies, radial distribution functions, and radii of gyration. Because such observables as the potential energy and radial distribution functions are diagonal in the position representation and therefore do not directly depend upon the fractional exponent, they can be measured as in the case. Appropriate care must simply be taken for cases because these averages are computed from distributions with heavy tails.

Computing quantities based upon derivatives of the partition function, such as the kinetic energy, however, is much less straightforward. As derived in Appendix III, expressions dependent on derivatives of the partition function are proportional to slightly more complex Fox-Wright distributions. The kinetic energy, for example, may be written as

where and denotes the average over links between beads. Note that the additional complexity results from the appearance of factors of as opposed to in the gamma functions. As a result of this additional complexity, the series is less manageable. It is therefore more convenient to evaluate the thermodynamic estimator of the average kinetic energy based upon the formula and to use a three-point forward difference formula in to estimate the derivative (see Section V).

## Iii Results and Discussion

### iii.1 Fractional Density Matrix Distributions

In this section, we compute the Fox-Wright density as described above to analyze its functional form and the performance of our sampling technique. In Figure 3, we illustrate the overall behavior of the probability density function that represents the uni-dimensional fractional density matrix. The density exhibits “heavy” algebraic tails for fractional order . Compared to the normal case, the fractional cases exhibit a larger probability of both small and large jumps, and a lower probability for medium jumps. In Figure 4, we look more closely at the tail behavior. As decreases, the related probability density functions exhibit heavier and heavier tails.

In Figures 3 and 4, we have also plotted the density for as a curiosity. Of course, -stable distributions are only defined for ; this illustration shows that for , takes negative values and is therefore not a density at all. In particular, while corresponds to superduffision, does not correspond in any way to subdiffusion. This may be surprising at first, but it has the following heuristic explanation. Unlike superdiffusion, one cannot obtain subdiffusion simply by modifying the spatial jump density of Brownian motion. In superdiffusion, the normal jumps are replaced by “faster”, infinite variance jumps. If one attempts to replace the jumps of Brownian motion by “slower jumps,” such a process would necessarily have finite variance and thus, by the central limit theorem, reduce to Brownian motion. Therefore, it is not surprising that for , is meaningless as a density. The only way to obtain subdiffusion is to introduce waiting times into the process, which results in a time fractional derivative.Meerschaert and Sikorskii (2011)

In Figure 5, we plot the same densities on a log-log scale to better understand their tail behavior. We can see from this that the density enters a power-law regime fairly early, for .

Finally, in Figure 6, we illustrate the accuracy of the discrete inverse transform sampling method described in Section II.5. This figure compares histograms representative of the probability density functions obtained for and obtained using our inversion algorithm. Note that our inversion method preserves the heavy tails that are central to these distributions.

The numerical costs of our inversion process are presented in Table 1. We see how expensive it is to evaluate the density : it takes about 1s to evaluate the density at 25 different points. On the other hand, for the grid used, it takes binary search 1s to generate 6000 samples. Moreover, since binary search is logarithmic in the number of discretization points, similar performance may be expected even for much finer grids. Therefore, if the CDF can be initialized to the desired accuracy, binary search is superior to Newton methods when a very large number of samples is required, i.e., orders of magnitude more than the number of discretization points. A Newton method would require several additional evaluations of for each sample; on a single core, this would generate samples per second, compared to the roughly thousands of samples per second for binary search.

### iii.2 Fractional Kinetic Energy

Using the forward difference scheme described in Appendix III, we analyze the functional form of the kinetic energy as a function of the interbead distance squared for a range of ’s. Interestingly, whereas the kinetic energy decreases linearly with interbead distance for , it first decreases, then increases, and finally plateaus for smaller , as depicted in Figure 7. Similar functional forms are observed for varying in Figure 8.

The plateaus observed in the kinetic energy may be understood by considering the power-law behavior of our Fox-Wright distribution for large argument. In general, the density matrix connecting two consecutive imaginary times, may be approximated as (see Appendix III)

 ρi(τ)=(21/αm1/α√πτ1/αℏα)3N× (36) [∞∑n=0Γ(1/α+2/αn)Γ(1/2+n)[−(→Ri−→Ri+1)222/αm2/α4τ2/αℏ2]n/n!].

Setting all parameters except , , and equal to 1 in the above, we can write

 ρi(τ)∼(m1α)3NFα(m1α√(→Ri−→Ri+1)2), (37)

where is the Fox-Wright distribution given by Equation (25). We note that Equation (25) has its argument squared in the Taylor series (see Equation (24)), so to write in terms of we must take the square-root of the expression

 [−(→Ri−→Ri+1)222/αm2/α4τ2/αℏ2] (38)

to use as an argument of . From Figure 5 above, we know that for large ,

 Fα(x)∼x−P,

for some that depends on . Thus,

 ρi(τ)∼(m3Nα)(m1α√(→Ri−→Ri+1)2)−P. (39)

As a result,

 dρi(τ)dm∼(m3Nα−1)(m1α√(→Ri−→Ri+1)2)−P − P(m3Nα)(m1α√(→Ri−→Ri+1)2)−P−1√(→Ri−→Ri+1)2(m1α−1)

and so

 mτρi(τ)dρi(τ)dm ∼ m−1τ−Pm1/α−1τ(m1α√(→Ri−→Ri+1)2)−1 (41) ×√(→Ri−→Ri+1)2 = m−1τ−Pm−1τ(m−1α)

which is a constant value. Thus, will initially decay and then approach this value for large .

As may be expected on a conceptual basis and is confirmed in the above plots, the average interbead distance for which the kinetic energy reaches its minimum shifts to larger values with increasing and . This is because larger imply more diffuse paths, while larger correspond to lower temperatures and therefore larger de Broglie wave lengths. The presence of a minimum in the kinetic energy is a strictly fractional feature which appears for infinitesimally smaller than two and suggests that fractional particles may exhibit intriguing behavior as they attempt to minimize their potential and kinetic energies – which have competing minima – at once.

### iii.3 Fractional Free Particle PIMC Results

Before presenting results for He, we first demonstrate our algorithm by simulating the free-particle Hamiltonian. Our interest in the free-particle Hamiltonian stems from the fact that it is amenable to analytic solutions and free-particle PIMC serves as the bedrock for interacting simulations. More specifically, we consider the Hamiltonian

 ^Hfree,α=−ℏα(2m)α/2∇α=−ℏαDα∇α (42)

for to avoid sampling from negative distributions. In order to make direct comparisons with our helium results presented in the next section, we set our mass equal to that of helium. In both sections, we moreover simulate particles at a density of bohr at all of the temperatures and values presented.

As a qualitative check on our simulations and underlying theory, we first consider the conformations of the polymers representing our particles. One of the principle motivations for this work was the fact that a fractional Laplacian for should lead to significantly more diffuse particles than in the typical case. As illustrated in Figure 4, this is because the fractional density matrix probability distributions possess heavy tails that should permit more frequent sampling of larger interbead distances, .

Figure 9 demonstrates that this is indeed the case: as decreases, the average square of the interbead distance, , for fixed and number of beads steadily increases. As may be expected, this increase is most dramatic at the low temperatures at which, in the absence of a potential, increased delocalization should be favored. Similar trends may be observed in the particles’ average radii of gyration, , as depicted in Figure 10.

Example paths illustrating these features for varying temperatures and fractional exponents are depicted in Figure 11.

What can be gleaned from these plots is that measures of both the average interbead distance and the radius of gyration are accompanied by increasingly large error bars as decreases. This can be expected as the variance of an -stable law is infinite for , with the mean itself becoming infinite for (we only consider the case). Although we have not implemented this in our current code, these variances may be made finite by sampling tempered distributions Meerschaert, Sabzikar, and Chen (2015) at the cost of approximating the distributions’ tails.Kleinert (2009)

Although related to the interbead distance, a more quantitative measure of our fractional paths is the average kinetic energy. In Figure 12, we plot the simulated PIMC kinetic energy vs. temperature for a range of fractional exponent values.

As may be expected from the increased interbead distances, the average kinetic energy increases with decreasing . The magnitudes of the kinetic energy may be rationalized by using the fractional version of the equipartition theorem: , which can be readily derived from Equation (III.2) in the high temperature limit. The larger kinetic energies for smaller suggests that there is a comparative kinetic energy penalty for fractional cases. While not illustrated here, if the average kinetic energy obtained from PIMC is plotted against the the average interbead distance squared, functional forms possessing minima and plateaus similar to those depicted in Figures 7 and 8 are obtained.

### iii.4 Fractional He-4 PIMC Results

Building upon our free particle simulations, we also simulate fractional He to demonstrate our PIMC algorithm in the presence of interactions. He is a particularly illuminating example because, even in the case, its particle paths become so diffuse at low temperatures that it can form a superfluid.Ceperley and Pollock (1986); Pollock and Ceperley (1984) In the following, we model He with the Aziz potential.Aziz et al. (1979) While the average interacting interbead distance remains similar to the non-interacting interbead distance, the average interacting radius of gyration seems to decrease significantly when compared to the non-interacting radius of gyration. This implies that the interacting particles are more localized (see Figure 14). This is also borne out in the decreased average kinetic energy depicted in Figure 15.

Irrespective of the increased localization due to the potential, the overarching trend that delocalization is favored with decreasing persists.

Interestingly, for He, we do not observe a pronounced competition between minimizing the kinetic and potential energies. This competition may be more pronounced, leading to more intriguing behavior, for interparticle potentials with deeper minima.

## Iv Summary and Outlook

In this work, we have presented a methodology that generalizes the Path Integral Monte Carlo algorithm to fractional Hamiltonians, thereby enabling the computational exploration of fractional Hamiltonians with potentials that are not readily amenable to analytical treatments. We have derived and shown how to sample fractional free particle density matrices by developing a novel approach to sampling Fox-Wright functions. This overall approach may be generalized beyond fractional path integrals to the sampling of density matrices of non-Gaussian forms. We have furthermore employed our algorithm to explore the impact of a fractional Laplacian on free-particle and He path observables, such as the radius of gyration and the average interbead distance. We have demonstrated that the fractional kinetic energy possesses a pronounced minimum as a function of average interbead distance for fractional exponents less than two () and that this minimum shifts to larger interbead distances with decreasing fractional exponents. As such, fractional particles become increasingly more diffuse with decreasing fractional exponents as a consequence of sampling a fractional density matrix distribution with heavier tails than the usual Gaussian distribution. Because the onset of phenomena related to condensation heavily depends on the diffusivity of particle paths, our findings suggest that fractional Hamiltonians may manifest intriguing superfluid and Wigner crystallization phase transitions. Such transitions have only previously been studied for the fractional Schrödinger Equation based upon mean field theories.Kleinert (2012) Fractional Path Integral Monte Carlo will enable a more complete understanding of such phenomena.

Before the effects of fractional operators on condensation phenomena may be explored, however, several complications revealed by our work must be resolved. First and foremost, because the conventional staging algorithmCeperley (1995); Sprik, Klein, and Chandler (1985) cannot readily be applied to non-Gaussian distributions, we sampled our partition functions using highly inefficient single bead moves in this work. While, for the potentials explored, we verified that our simulations were able to converge to their equilibrium distributions, we would not expect this to be the case for more general potentials. It is furthermore easy to imagine that single bead moves would be suboptimal for sampling permutations among particles with quantum statistics. Future work therefore necessitates the development of generalized multi-bead moves. Fourier path sampling may permit the generality we require.Doll, Coalson, and Freeman (1985); Freeman, Coalson, and Doll (1986); Krauth (2006)

As Figures 3 and 4 illustrate, the density matrix distributions sampled for in our algorithm are furthermore significantly more complex than the Gaussian distribution sampled conventionally. For , these distributions possess heavy tails, while for , these distributions can be negative. The discrete inverse-transform sampling method we presented can accurately sample the heavy tails, but there remain theoretical questions about the effect of heavy tails in the theory. Namely, whether the infinite variance (for ) and infinite mean (for ) of the fractional density matrix requires the use of other notions of central tendency (such as median) for observables. This is a fascinating question for future study. In this direction, it may be fruitful to explore “tempered-fractional quantum mechanics,” which has not been studied to our knowledge. This would conceivably employ tempered Levy distributions instead of pure heavy-tailed distributions, and temperated fractional derivatives instead of standard fractional derivatives.Meerschaert, Sabzikar, and Chen (2015) Tempered distributions enjoy the heavy tail to a certain argument, but then decay exponentially to avoid issues with infinite variance.

It may also be of interest to consider the physics of the time-fractional Schrodinger Equation,Naber (2004) and even the space-and-time-fractional Schrodinger Equation.Wang and Xu (2007) These equations introduce waiting times in addition to long jumps. However, unlike the space-fractional equation, the properties of the time-fractional Schrodinger equation suggest that it cannot be used as a probabalistic theory of Quantum Mechanics.Dong and Xu (2008b); Laskin (2017) Thus, more work is needed before any implementation of a time-fractional model is justified.

In summary, we have proposed a novel fractional Path Integral Monte Carlo algorithm. Our algorithm provides a clear numerical path toward unraveling the fractional physics of interacting systems, much of which could only be approximated analytically in the past. We look forward to employing this algorithm to explore where our fully interacting predictions differ from previous approximate analytical solutions and to the study of such novel phenomena as fractional superfluidity.

###### Acknowledgements.
The authors would like to dedicate this work to the late Jimmie Doll. Jimmie initially convened the authors to brainstorm about this topic, but unfortunately, did not live to see this work to its fruition. The authors would also like to thank George Karniadakis, who suggested to M.G. that he explore this topic and contact Jimmie. Finally, we thank Paul Dupuis, Matthew Harrison, and Mark Meerschaert for valuable discussions. M. G. would like to acknowledge support from the OSD/ARO/MURI grant “Fractional PDEs for Conservation Laws and Beyond: Theory, Numerics and Applications (W911NF-15-1-0562).” Part of this research was conducted using computational resources and services at the Center for Computation and Visualization, Brown University.

## V Appendix I: Separability of the Fractional Density Matrix in Multiple Dimensions

In this work, we have derived and demonstrated how to sample the one-dimensional form of the free fractional density matrix. Here, we reaffirm that, just as in the case, a multi-dimensional free density matrix is simply a product of one-dimensional free density matrices.

Recalling Equation (5), the free multi-dimensional density matrix may be written as

 ρ(→R,→R′;τ) = ⟨→R|e−τ^K|→R′⟩ (43) =</