Although disorder (especially strong like substance amorphization) in a solid is usually considered to be a trouble for its possible device applications, its constructive properties have become increasingly clear in recent years as they give an additional possibility to fine tune their physical properties. Specifically, by varying the type and concentration of different imperfections (like point or extended defects), we can customize the characteristics of such material to meet specific requirements, needed, for instance, in electronics. Disordered semiconductors and dielectrics like organometallic halide perovskites^1,2 are widely used in photovoltaic cells, light-emitting devices, and nanolasers^3,4. The above disorder influences phonon and electron spectra of a substance, leading to the distribution of the internal magnetic, electric and elastic fields. This is the case for so-called multiferroics, where ferroelectric and magnetic orders coexist⁵. It is well-known (see, e.g.^6,7) that strong disorder is well described by Non-Gaussian probability distributions, having so-called long tails. In other words, these distributions decay at infinities (sometimes much) slower, that Gaussian one. The most prominent example here is so-called Lévy flights^,,,,,13, which are Markovian random processes whose probability density functions (pdf) are Lévy stable laws, characterized by the Lévy index \(0<\mu \le 2\). In the infinite space, it is profitable to define the corresponding pdf in terms of its characteristic function, i.e. Fourier transform. The pdf of a typical Lévy flight is usually determined by the fractional Fokker–Planck (FP) Eq. (⁸). In dimensionless units (diffusion coefficient and particle mass are set to unity) it reads

where \(|\Delta |^{\mu /2}\) is an n-dimensional fractional Laplacian (Riesz fractional derivative^6,7)

which at \(\mu =2\) yields the ordinary one^6,7. Here \(\Gamma (x)\) is \(\Gamma \)-function¹⁴, \(\nabla \) is (also n- dimensional) gradient operator and \(U(\mathbf{r})\) is an external potential, in which the anomalous diffusion occurs¹⁵. The operator (Eq. 2) is spatially nonlocal, which is a source of memory effects in strongly disordered systems.

It had been shown by Laskin, that the substitution of the Gaussian measure by the Lévy one in the Feynman path integrals generates so-called fractional quantum mechanics^16,17, where an ordinary Laplacian in the Schrödinger equation is substituted with fractional one (Eq. 2). In context of the above disordered solids, it is important to consider the influence of disorder onto the spectra of their collective atomic vibrations, which constitute well-known phonons in an ordered substance. As the latter vibrations are described by the quantum harmonic oscillator model, it is reasonable to consider the fractional analog of real 3D quantum oscillator problem.

In the present paper we shall find the spectrum of 3D fractional oscillator for arbitrary index \(0<\mu \le 2\). As this problem resides on the whole 3D space, we shall solve it using Fourier transformation to the momentum space. In latter space the problem becomes local as the Fourier image of the operator (Eq. 2) is simply \(k^\mu \) (\(k\equiv |\mathbf{k}|\), where \(\mathbf{k}\) is the momentum vector), while the oscillator potential gives the ordinary Laplacian in the momentum space. In other words, the initial fractional Schrödinger equation gives ordinary one in \(\mathbf{k}\) space. The only difference is that the potential in this new equation will be \(k^\mu \). This implies that to solve our problem in momentum space, we can use well developed numerical and variational techniques.

In the units, where \(\hbar =m=\omega =1\) (m and \(\omega \) are the mass and frequency of oscillating particle respectively), the spectral problem for 3D fractional oscillator reads

where the operator \(|\Delta |^{\mu /2}\) is defined by the expression (2) and \(\psi _{nlm\mu }(x)\) is the eigenfunction of a fractional quantum harmonic oscillator having the eigenenergy \(E_{nl\mu }\) for any specific \(\mu \) value. In general 3D case, as it is typical for the quantum-mechanical problems in central force potentials, the eigenfunction \(\psi \) and eigenenergy E of the problem (4) depend on three quantum numbers n, l, m which are principal, orbital and magnetic quantum numbers respectively^18,19.

As we discussed above, the most profitable method to solve this problem is to pass to the momentum space, which permits to reduce the integral equation (4) to the differential one. Namely, as it is well-known that the combinations \(x^2\psi (x,y,z)\), \(y^2\psi (x,y,z)\) and \(z^2\psi (x,y,z)\) (here we suppress the indices \(nlm\mu \) for brevity) in momentum space render as second derivatives \(\frac{\partial ^2}{\partial k_x^2}\psi (k_x,k_y,k_z)\), \(\frac{\partial ^2}{\partial k_y^2}\psi (k_x,k_y,k_z)\) and \(\frac{\partial ^2}{\partial k_z^2}\psi (k_x,k_y,k_z)\), the potential term in the equation (4) renders simply as an ordinary Laplacian \(\Delta _\mathbf{k}\equiv \frac{\partial ^2}{\partial k_x^2}+\frac{\partial ^2}{\partial k_y^2}+\frac{\partial ^2}{\partial k_z^2}\) in momentum space. At the same time, the fractional Laplacian in momentum space becomes simply \(k^\mu \), where \(k\equiv =\sqrt{k_x^2+k_y^2+k_z^2}\) is a vector \(\mathbf{k}\) modulus. Latter definition immediately implies that \(0\le k< \infty \). This means that in momentum space, the equation (4) renders as

where \({\mathcal H}_\mathbf{k}\) is the system Hamiltonian in momentum space. One note is in place here. Namely, the equation (5) shows that at \(\mu <2\), the famous orbital 3D oscillator degeneracy (see Ref.¹⁸ for details) is lifted so that the eigenenergy E starts to depend on the orbital quantum number l. At the same time, as in the fractional case of \(\mu <2\), the time-reversal symmetry remains (it is lifted, for instance, in the external magnetic field¹⁸), we still have degeneracy with respect to the quantum number m. Namely, the spectrum is invariant as \(m \leftrightarrow -m\). That is why, in equation (5), we put index m explicitly in the wave functions, which are different for different m, although the energy E does not depend on m.

As for any central potential the wave functions are invariant with respect to rotation around the origin^18,19, this is also the case for our 3D fractional oscillator. This implies that the wave functions \(\psi _{nlm\mu }(\mathbf{k})\) admit the separation of angular and radial (in \(\mathbf{k}\) space) variables in the form

where \(Y_{lm}(\theta , \varphi )=A_{lm}P_l^{|m|}(\cos \theta )e^{im\varphi }\) are spherical harmonics^14,18. Here \(P_l^{|m|}(\cos \theta )\) are associate Legendre polynomials of the first kind and \(A_{lm}\) is the corresponding normalization coefficient, which is obtained from the following normalization condition^14,18

where \(d\Omega =\sin \theta d\theta d\varphi \).

It is well-known from “ordinary” quantum mechanics^18,19 that 3D quantum oscillator problem for \(\mu =2\) admits the factorization of the wave function \(\psi (x,y,z)=\psi _x(x)\psi _y(y)\psi _z(z)\) (here we also suppress the lower indices for a moment), which permits to reduce the problem to the three equations for 1D oscillator, containing three different energies, say \(E_x\), \(E_y\) and \(E_z\). In this case, the initial energy \(E=E_x+E_y+E_z\). After problem being in such a way factorized, the eigenenergy can be derived in the well-known form^18,19

where \(n=n_x+n_y+n_z\) with \(n_{x,y,z}\) being the quantum numbers, corresponding to \(E_{x,y,z}\) respectively. Such factorization in the “ordinary” case \(\mu =2\) reduces the problem of the 3D oscillator to the combinatorial problem with respect to its 1D “components”. In other words, the quantum numbers nlm of the resulting 3D problem comprise different combinations (obeying certain combinatorial rules, see Refs.^18,19 for details) of the initial \(n_x\), \(n_y\) and \(n_z\) numbers.

It had been shown in¹⁹, that the possibility of such factorization at \(\mu =2\) is based on additivity of both Laplacian operator (kinetic energy) and the potential \(\sim (x^2+y^2+z^2)\). At the same time, in \(\mathbf{r}\) space (Eq.4), the fractional Laplacian (Eq. 2) is not additive. That is to say, it does not consist of a sum of the terms, containing separately x, y and z variables. This fact is “transferred” to \(\mathbf{k}\) space, where the potential \(k^\mu \equiv (k_x^2+k_y^2+k_z^2)^{\mu /2}\) is factorable for \(\mu =2\) (ordinary case) only.

This signifies the difference between fractional \(\mu <2\) and ordinary \(\mu =2\) 3D oscillator problems. This difference implies that for our fractional case we should proceed directly in 3D case. This, however, does not give an essential complication as problem still reduces to 1D one. Namely, after the decomposition (Eq. 6), we obtain following fractional 1D Scrödinger equation for the radial part \(R_{nlm\mu }(k)\)

Here we also suppress the indices \(nlm\mu \) in the functions R and energy E for a moment. In the quantum mechanics of the central force systems, it is also customary (see, e.g.^18,19) to pass to the function

for which the Eq. (9) assumes particularly simple form

which resembles that for 1D system. This form will be used below to construct trial functions for variational method as well as for numerical calculations, where the form (Eq. 11) permits to avoid divergencies at \(k=0\). We note here, that if for the states with \(l=0\), the function R(k) has a maximum at \(k=0\), the corresponding function \(\chi \) is zero by obvious reason.

As we mentioned above, the Eq. (11) can be readily solved numerically. On the other hand, even approximate analytical solution like variational one is much more desired as it can represent the useful orthogonal set, which can be further employed for the solutions of the problems, not obligatory dealing with fractional quantum mechanics. For instance, it could be well used to look for the solutions of the fractional FP equation (1) in the from of the expansion over the infinite sets of orthogonal functions. That is why we will solve the radial equation (11) variationally and confirm our results by direct numerical solution.

Potential in the Schrödinger equation (9) at \(l=0\) (a) and \(l=1\) (b). Lévy indices \(\mu \) are coded by colors (legend in (a)). Mind the different vertical scales in (a,b).

Full size image

The Eqs. (9), (11) represent Schrödinger equation with the effective 1D potential \(|k|^\mu +l(l+1)/k^2\). The plots of this potential at different \(\mu \) and l are shown in Fig.1. It is seen from Fig.1a that for \(l=0\) the problem is identical to that for 1D fractional oscillator²⁰, while for \(l \ne 0\) (Fig.1b), the centrifugal force drives it outwards. The only difference with the case of 1D oscillator is that now the problem is defined on positive semi-axis \(0<k<\infty \) as k is now the modulus of a momentum vector. As it is seen from Fig.1b, the centrifugal term in the potential \(l(l+1)/k^2\) is dominant at \(k \rightarrow 0\), while at \(k \rightarrow \infty \) the main contribution comes from \(k^\mu \). This means that the whole spectrum of 3D fractional oscillator remains to be discrete similar to the “ordinary” case \(\mu =2\). Note that at small \(\mu =0.2\) the potential begins to grow at very large k, while at \(\mu =0\) the potential goes to 1 at \(k \rightarrow \infty \). The same situation at \(l=0\) gives constant value one for \(\mu =0\), see Fig. 1a. Our numerical analysis shows that in 3D case despite the contribution of centrifugal term, at \(\mu =0\), the whole oscillator spectrum collapses into one value, see below. This property is similar to that in 1D case²⁰. Note that at \(l=0\) and \(\mu <1\) the potential at \(k \rightarrow 0\) the potential has nonanalytical behavior with infinite first derivative. At \(l \ne 0\) this asymptotics is “killed” by the centrifugal term.

It is well-known that the variational principle of quantum mechanics¹⁸ permits to find the approximate analytical expressions for the wave functions and eigenenergies of corresponding Schrödinger equations. This is especially true for the states of discrete spectrum, which are spatially localized. As usually, the “correct” (i.e. that which minimizes corresponding energy functional with minimal number of variational parameters) set of trial functions should be obtained on the class of functions, having the asymptotics, dictated by the corresponding differential equation.

Note that for the above variational principle to work, it is sufficient even to use the trial functions obeying the correct boundary conditions rather than have the “prescribed” asymptotics. In our case, this means that the trial functions have to be zero at \(k \rightarrow \infty \). Such “simple” trial functions can be constructed on the base of the Gaussian one, inherent in an “ordinary” (i.e. that at \(\mu =2\)) 3D oscillator. In this case, the whole variational spectrum can be thought of in the form of “Hermite polynomials with unknown coefficients”, where the latter coefficients are to be found from the orthogonality conditions, see below. In the next section, we shall compare the relative errors of the variational approaches, based on “fractional” (i.e. with predefined asymptotics) and “Gaussian” trial functions.

To construct the “fractional” trial functions, we should analyze the large k asymptotics of the \(\chi (k)\), given by the radial (in \(\mathbf{k}\) space) Schrödinger equation (11). For that we observe that at large k we can neglect all the terms in the square brackets in Eq. (11) except \(k^\mu \). This yields the equation for large k in the form \(d^2\chi (k)/dk^2=k^\mu \chi \). The spatially decaying solution to this equation is proportional to \(k^{1/2}K_\nu (u)\), where \(\nu =1/(\mu +2)\) and \(u=2\sqrt{2}|k|^{1+\mu /2}/(\mu +2)\)²¹. Here \(K_\nu (x)\) is so-called MacDonald function with the following large x asymptotics \(K_\nu (x \rightarrow \infty )\approx (\pi /(2x))^{1/2}e^{-x}\)¹⁴. Substitution of latter asymptotics into the above expression generates following wave function asymptotics at all admissible \(\mu \)’s

Here we once more suppressed indices nlm. In the case of \(\mu =2\) (ordinary quantum oscillator) the asymptotics (Eq. 12) reproduces well known result for “ordinary” 3D quantum oscillator¹⁸. It is seen that for all \(0<\mu \le 2\) the asymptotics (Eq. 12) decays sufficiently fast for corresponding integrals to be convergent. Namely, in worst case of \(k=0\) \(\psi (k \rightarrow \infty ) \sim e^{-k\sqrt{2}}\), i.e. decays exponentially. Good localization of the wave functions in k space yields their localization in coordinate space by Riemann–Lebesgue lemma, see e.g.²². Our analysis show that at small k the functions R(k) have the same asymptotics \(R \sim k^l\) (and accordingly \(\chi \sim k^{l+1}\)) as in ordinary case, corresponding to \(\mu =2\)¹⁸.

To classify the states of the 3D fractional oscillator according to the quantum numbers nlm, we notice that although the Eq. (9) resembles that for the ordinary 3D oscillator, it is not equivalent to it as we have the term \(k^\mu \) instead of \(k^2\) in the potential. The former term “spoils” the neat classification^18,19 of the states of 3D quantum mechanical oscillator so that in the fractional case \(\mu <2\) there is no obvious way to introduce the principal quantum number n. We recollect that for ordinary 3D oscillator \(n=2n_r+l\) (see, e.g.¹⁹), where \(n_r\) is so-called radial quantum number, defining the number of nodes of the radial function for a given l. That is why to classify the states of the oscillator in our 3D fractional case we begin with general consideration of a quantum particle, moving in a central force field. Namely, to apply the oscillation theorem to the radial part \(\chi (k)\) in Eq. (11), we arrange the eigenenergies corresponding to a given l in ascending order. It had been shown^18,19, that the lowest eigenenergy for a given l corresponds to a nodeless function so that we assign \(n_r=0\) to this state. The next (excited) state for the same l has one node, i.e. \(n_r=1\) etc. This means that we may classify the states of our fractional oscillator by essentially two quantum numbers - \(n_r\) and l. To be specific, for \(l=0\) (s - state; here we use the common notations of the states with the given l, see, e.g.¹⁸) we have the whole sequence \(n_r=0,1,2,...\); the same is true for \(l=1,2,3,..\). Introducing the principal quantum number \(n=n_r+1\), we obtain following explicit classification of the states of the fractional quantum oscillator: 1s (\(n_r=0\), hence \(n=1\), \(l=0\)), 1p (\(n=1\), \(l=1\)), 1d (\(n=1\), \(l=2\)),... 2s, 2p, 2d,... etc. In this case, the commonly known degeneracy of the states 1d and 2s, 1f and 2p etc, inherent in “ordinary” 3D oscillator, is lifted so that all states with different \(n_r\) and l have different energies. This classification will be used below to construct the trial functions.

To construct the above functions, we need the asymptotics at small and large k as well as the information about the number of nodes, i.e. the numbers \(n_r\) and l. Specifically, in terms of functions \(\chi (k)\), we have

Here, \(A_{nl}\) are normalization constants, found from the “1D-like” conditions

Also, \(a_{nl}\) are variational parameters, while the constants \(a_i\), \(b_i\) and \(c_i\) (here \(i=1-4\)) are found from the orthogonality conditions

where \(\delta _{jj'}\) is Cronecker delta. We note that the trial functions \(\chi _{nl}\) for different l (for instance \(\chi _{1p}\) and \(\chi _{2d}\)) are orthogonal automatically because of orthogonality of spherical harmonics^14,18,19. That’s why the integration in Eqs. (16) and (17) runs over k only. This implies, that only the functions with different \(n's\) (and the same l) should be orthogonalized. Namely, in the particular case of the above trial functions we have explicitly

The orthogonalization procedures Eqs. (17), (18) along with normalization of the functions \(\chi _{nl}\) determine uniquely the coefficients \(a_i\) through variational parameters \(a_{nl}\).

We have explicitly

Here \(\Gamma (x)\) is the \(\Gamma \)-function¹⁴. The other normalization coefficients can be calculated, but their explicit expressions become progressively more cumbersome.

Orthogonality condition \(<\chi _{1s}\chi _{2s}>=0\) yields

Substitution of the expression (20) into the normalization coefficient \(A_{2s}\) Eq. (19b) generates following explicit expression for the function \(\chi _{2s}\)

where \(\lambda \) is determined by the expression (20) and

Now, after variational determination of the parameter \(a_{1s}\) (see below), the parameter \(\lambda \) (20) starts to depend on the single parameter \(a_{2s}\), which can also be determined variationally. Having determined the parameters \(a_{1s}\) and \(a_{2s}\), we then use the orthogonality condition \(<\chi _{1s}\chi _{3s}>=0\) to express both \(c_2\) and \(b_2\) through \(a_2\). After substitution of the above found \(a_{1s}\) and \(a_{2s}\) into the corresponding normalization coefficient, the function \(\chi _{3s}\) starts to depend on \(a_{3s}\) only. Such “chain rule” can be applied to determine uniquely all higher order functions, corresponding to \(l=0\) (s-state), \(l=1\) (p-state) and for any arbitrary n and l.

The expressions (13)–(15) show the algorithm of construction of trial functions. This algorithm is essentially the same as that for “ordinary” case of \(\mu =2\). Namely, the arbitrary trial function \(\chi _{nl}\) is the product of the factors \(k^{l+1}\) (small k asymptotics), \(e^{-a_{nl}k^{1+\mu /2}}\) (large k asymptotics) and a polynomial of degree \(n_r=n-1\) with unknown coefficients. These coefficients, in turn, are determined from the orthogonality conditions (17) by the above “chain rule”.

As usually, the variational solution of the Schrödinger equations should minimize the energies

where \(\mathcal {H}_\mathbf{k}\) is the Hamiltonian (5) and \(W_{nl\mu }\) are the variational approximations of the eigenenergies \(E_{nl\mu }\). Normally \(W_{nl\mu } \ge E_{nl\mu }\). As the radial parts \(R_{nlm\mu }(k)\) of the wave functions (6) are well localized (see asymptotics (12)), the expression (23) could be rendered (by the integration by parts) to the more convenient form

where \(R'=dR/dk\). After substitution of the trial functions (13), (14) and (15) into the integral (24), with respect to the above “chain rule”, we find the variational parameters \(a_{nl}\) from its minimization. Such a procedure can be done for all wave functions of higher excited states, giving the approximate spectrum of the operator (5). Note that the substitution of found \(a_{1s}\) into \(W_{1s\mu }\) (24) gives the approximate value of the ground state energy for all \(\mu \), the same with, say, \(a_{2p}\) gives the energy of the excited state \(W_{2p\mu }\) and so on for higher eigenenergies. Although the expressions for higher excited states become extremely cumbersome (see above), the variational method permits to obtain the approximate analytical expressions for the eigenvalues and eigenfunctions of the operator (Eq. 5) for all admissible \(0<\mu \le 2\).

Substitution of the trial function \(\chi _{1s}\) into the integral (Eq. 24) with subsequent minimization over \(a_{1s}\) yields

which is similar to the expression for 1D fractional oscillator²⁰. Further substitution of this value to the expression (24) generates the approximate value of the ground state energy for arbitrary \(\mu \)

It is seen that \((W_{1s\mu })_{min}\) gives correct value 3/2 of the ground state energy for \(\mu =2\), corresponding to the ordinary quantum oscillator with the spectrum \(E_n=n+3/2\) in our units. Below we shall see, that for the case \(\mu =0\) all the spectrum shrinks into a single value \(E_0=1/2\), which is also obtained correctly from the expression (26). The dependence (Eq. 26) will be plotted below and compared with the numerical solution.

The same procedure with \(\chi _{1p\mu }\) gives that \((a_{1p\mu })_{min}=(a_{1s\mu })_{min}\), which is given by Ex. (25). The variational expression for the energy of the 1p state reads

The expression (27) shows that at \(\mu =2\) we obtain the correct answer 5/2. We recollect here that the energy of “ordinary” 3D quantum oscillator is determined by the expression (8), where \(n=2n_r+l\), see above. As for 1p state we have \(n_r=0\) (no nodes) and \(l=1\) so that \(n=1\) and the energy \(E_{1p}(\mu =2)=1+3/2=5/2\). At the same time, at \(\mu \rightarrow 0\) we have removable divergence \(\lim _{\mu \rightarrow 0}(2\mu ^{-\mu /2})=\lim _{\mu \rightarrow 0}\exp (-\mu \ln (2\mu ))=1\), while the rest of the expression (27) gives \(E_{1p,\mu =0}=6\Gamma (3)/\Gamma (5)=1/2\), i.e. once more the correct answer. The dependence (27) will also be plotted below and compared with numerical results. We will see that the expressions (26) and (27) give very good approximate expressions for corresponding energies of a 3D fractional quantum oscillator for all admissible \(\mu \). Within the suggested variational approach, the energy of any desired level \(E_{nl\mu }\) can be evaluated analytically, although the calculations for higher excited states become very cumbersome, see above.

(a) Reports the lowest state eigenenergies (indicated near the curves) as functions of Lévy index \(\mu \), calculated numerically (solid lines) and variationally (dashed lines, “fractional” functions (13)–(15)). (b) Displays wave functions of the lowest states, shown in the legend. The numerical and variational plots in (b) are similar in the scale of the plot. (c) The 1s state energies, calculated variationally (expressions (26) and (31)) as well as numerically. The curves are indistinguishable in the scale of the plot. Inset in (c) reports the “Gaussian” - numerical (black curve) and “fractional”—numerical (red curve) relative errors in the entire \(\mu \) domain. (d) Portrays the comparison between variational (both types of the functions, shown in the legend) and numerical 1s wave functions for the same \(\mu \) values, as in (b). Insets report the relative errors for \(\mu =1.4\) and 0.4 as the wave function for \(\mu =2\) is exact. Black curves correspond to “Gaussian”—numerical and red ones to “fractional”—numerical relative errors as the functions of momentum k. Figures near curves in (b,d) correspond to Lévy indices \(\mu \).

Full size image

The easiest way to solve our problem numerically is to use the Eq. (11) for the function \(\chi (k)\), which has only one divergency (in the centrifugal term \(l(l+1)/k^2)\)) at \(k=0\). The numerical solution of the Eq. (11) should be augmented by the boundary conditions \(\chi (0)=\chi (\infty )=0\). As usually for boundary value problems, our Eq. (11) can be solved numerically by its reduction to the eigenproblem for the corresponding finite-dimensional matrix. In this case, the energies E are the eigenvalues of the above matrix and the eigenvectors of the latter comprise the wave functions of the problem in momentum space. The transition to the usual coordinate space is accomplished by 3D inverse Fourier transformation. We note that the satisfactory accuracy of the numerical solution is achieved for typical matrix dimensions 10000 \(\times \) 10000, which makes the task quite computer intensive.

As we can reproduce numerically the spectrum of the problem for arbitrary \(\mu \), we are now in a position to compare the results of variational and numerical treatments. We begin with comparison of different classes of trial functions (see above), namely aforementioned “fractional” (13)–(15) and “Gaussian”, which we are going to construct below. As the consideration of a ground (1s in our case) state wave functions give the lower limit of the relative errors, here we consider just this function for “Gaussian” case. In terms of \(\chi \) functions (10), we have

where “G” in lower index stands for “Gaussian”, A is normalization coefficient and \(\sigma \) is variational parameter. We note that at \(\mu =2\) the function \(\chi _{G1s}\) is equivalent to \(\chi _{1s}\) (see (13)) with \(a_{1s}=1/(2\sigma ^2)\).

The normalization of (28) with the help of condition (16) yields

Substitution of Eq. 28 with respect to Eq. 29 into the expression for energy (Eq. 24) with its subsequent minimization over \(\sigma \) generates following minimizing wave function width

Further substitution of the ground state wave function (28), with respect to Eqs. 29 and 30, into the energy functional (24), yields the Gaussian analog of the variational ground state energy (Eq. 26)

where \(\sigma _{min}\) is determined by Eq. (30). At \(\mu =2\), the energy \((W_{G1s\mu })_{min}\) also gives correct value 3/2 of the oscillator ground state energy. The same is valid for \(\mu =0\), where \((W_{G1s\mu =0})_{min} =\Gamma (3/2)/\sqrt{\pi }=1/2\). Here we also have a removable divergency of the type \(\lim _{\mu \rightarrow 0}(\mu ^{-\mu /2})=\lim _{\mu \rightarrow 0}\exp (-(\mu /2)\ln \mu )=1\) in the first term of Eq. (31).

Having both “Gaussian” (31) and “fractional” (26) and (27) expressions for variational energy levels of our 3D oscillator, we are now in a position to compare them with corresponding numerical values. Such comparison is reported in Fig. 2a,c. As is the case for the variational method, in both panels the variational curves lie higher than numerical ones as variational energy should be larger than its exact (in our case numerical) value¹⁸. It is seen that the agreement is a little worse for the 2s state, which has one node. For the nodeless states with \(n=1\), the variational and numerical curves are indistinguishable in the scale of the plots. The quantitative error analysis for 1s state is reported in Fig. 2c. It shows that as both “Gaussian” and “Fractional” trial functions give exact results at \(\mu =0\) and \(\mu =2\), the relative error should be maximal somewhere inside this interval. The inset in Fig. 2c shows that the maximal relative error occurs at \(\mu \approx 0.5\) and does not exceed 0.6%. Latter error occurs for “Gaussian” trial function. This is because the latter function does not take into account the actual asymptotics (Eq. 12) of the wave functions.

The selected numerical wave functions in \(\mathbf{k}\)-space are reported in Fig. 2b,d for different values of \(\mu \). It can be checked that functions are normalized, i.e. they obey condition (16). Note that corresponding “fractional” wave functions (13) are almost indistinguishable from numerical ones in the scale of the plot in Fig. 2b. At the same time, in the Fig. 2d, we choose the vertical scale so that the (minute) differences between “fractional” (variational) and numerical 1s wave functions become visible. Panel (b) of Fig. 2 shows that functions for different \(\mu \) have different decay rates, which are dictated by the asymptotics (Eq. 12). For instance, the 2s function for \(\mu =0.4\) at \(k=5\) is not yet achieved its final part (compare to that for \(\mu =1.4\)), which occurs at \(k \approx 15\). For smaller \(\mu \), the spatial extension of the wave functions in \(\mathbf{k}\) space increases drastically. The functions for \(\mu \rightarrow 0\) in \(\mathbf{k}\) space are almost delocalized, which means that in coordinate space they are Dirac \(\delta \)-function-like. This behavior is qualitatively similar to that of 1D fractional quantum oscillator²⁰ as well as for 2D^23,24 and 3D²⁵ fractional hydrogenic problems.

The same tendency is seen in Fig. 2d. Here, the “Gaussian” trial function has stronger deviation from the numerical curve, while the “fractional” one goes closer to it. This difference is more visible at \(\mu =0.4\), which is close to the above value \(\mu \approx 0.5\), where the relative error in the energy is maximal. As numerical and variational curves are identical for \(\mu =2\) (the corresponding trial functions are exact), the relative error is reported for \(\mu =1.4\) and 0.4 only. The behavior of relative error reflects the well-known fact that different classes of trial functions can give a very good accuracy for the energy, while the approximation of the numerical wave function can be much worse. We note here that the values of relative error at both \(k \rightarrow 0\) and \(k \rightarrow \infty \) are not that reliable as in this case we are dealing with the difference of small (in the sense of machine precision) numbers, divided by those (small) numbers. That’s why we limit ourselves by \(k \approx 3\) as for \(k>3\) because of latter effect, the error grows to infinity. It is seen from the insets in Fig. 2d, that in the intermediate region of k values, the maximal relative error of the “Gaussian” trial function is around 10%, while that for “fractional” one is less than 1%. Our analysis shows that for the “fractional” trial functions with nodes (including 2s and higher) the average error is 2–3%, the larger error around 3% occurs for the higher excited states with many nodes. This shows that the variational expressions like (13)–(15) (with respect, of course, to (19a)–(19d)) can be regarded as pretty good wave functions in \(\mathbf{k}\) space. Such good variational approximation by the functions, having only one parameter may be because the energies are monotonous functions of \(\mu \). That said, we can safely assert that (in the intermediate range of momenta k, where the error can be calculated reliably, see above) with the accuracy not higher than 3% the analytical expressions (13)–(15) approximate the corresponding numerical curves.

Wave functions of 1s (a) and 2s (b) states in the coordinate space. The values of \(\mu \) are encoded by colors (legend in (a)).

Full size image

One more observation is in place here. Namely, the eigenenergies in fractional case \(\mu <2\) are smaller than those for “ordinary” case \(\mu =2\). This means that in the system, where the substitution of conventional Laplacian by the fractional one in Schrödinger equation is admissible (like in strongly disordered solids like amorphous ones, see, e.g.²³ for discussion), it is energetically favorable for an oscillator to “become fractional”.

As we have discussed above, the wave functions in the coordinate space can be obtained by the inverse Fourier transform. They are reported in the Fig. 3 for 1s (panel (a)) and 2s (panel (b)) states. It is seen from the Figure, that the curves are qualitatively similar to those in the momentum space. The regularities here are similar to those for 1D case²⁰, i.e. the slowest decay have the functions for \(\mu \rightarrow 2\) (“ordinary” quantum oscillator), while the fastest (like Dirac \(\delta \), see above) have the functions for \(\mu \rightarrow 0\) (\(\mu =0.4\) in Fig. 3), corresponding physically to the strongest disorder. This may be one more signature of the famous Anderson localization²⁶ phenomenon.

The obtained spectrum of the fractional quantum 3D oscillator can be used for the calculation of arbitrary observable characteristic of real physical systems. As an illustration, here we consider the so-called oscillator strength or dipole matrix element (the matrix element of a dipole moment operator \(\mathbf{d}=e\mathbf{r}\)) for the transition between the energy levels with different orbital quantum numbers l and \(l'=l \pm 1\). The latter relation stems from the selection rule, which gives zero (because of angular integration) for any combination of orbital quantum numbers, other than the latter one. The oscillator strength is important for the evaluation of such experimentally observable characteristics like the intensities of electric dipole transitions, which are the dominant effects of the interaction of a charge carrier in solids (electrons or holes) with the external electromagnetic field. One of the most important matrix elements of this type defines a transition between the states with \(l=1\) and \(l=0\). As here we are working in momentum space, the formula for the oscillator strength contains the matrix element \(\langle a|\mathbf{k}|b \rangle \) for the transition between states a and b. In our dimensionless units, the expression for the oscillator strength reads (see, e.g.^19,27)

where the transition occurs between states a (like 1p) and b (like 1s) with the eigenenergies \(E_a\) and \(E_b\) respectively. For the transition between the states 1p and 1s, the explicit expressions for matrix elements assume the form

where

Substitution of the matrix elements (Eqs. 33a–33c) with respect to numerically calculated \(R_{00}\) and \(R_{11}\) wave functions into the expression (32), generates the Lévy index dependence of oscillator strength components \(f_i\equiv f_{i, 1s-1p}\) (\(i=x,y,z\)), which is reported in Fig. 4. It is seen that the oscillator strength goes to zero in the case of strongest disorder, corresponding to \(\mu =0\), see Ref.²³ for details. We observe (see above) that as \(\mu \rightarrow 0\), both \(E_{1p}\) and \(E_{1s}\) go to the same value (as, actually, the whole spectrum) 1/2. This means that for \(\mu =0\), the denominator of the expression (32) becomes zero. Despite that, the whole expression (32) tends to zero because of much faster decay of dipole matrix elements (33a)–(33c) in this case. The origin of this stems from coordinate space, where the wave functions of all states look very much like Dirac \(\delta \)-functions. This means that the probability of optical transitions is much higher in relatively ordered regions, where \(\mu >1\). This fact may be useful in the theoretical description of the experimental data either in polymers, where chains can vibrate, or in disordered solids. Namely, the developed formalism can be well applied to the calculations of the properties of real physical systems, where disorder (like lattice imperfections and/or impurities) influences phonon and electron spectra of a substance, leading to non-Gaussian distribution of the internal electric, magnetic and elastic fields. A compelling example is the properties of the above multiferroics⁵. The non-Gaussian statistics due to disorder and frustration plays an important role in these substances^,,30 and it is challenging to apply the above formalism to explain their disorder (in phonon and magnon subsystems) phenomenologically in terms of introduction fractional derivatives to the corresponding Schrödinger equations. The dynamics of such “disordered phonons and magnons” can be described in terms of the above FP equation. Here, essentially the same variational approach can be used as there exist a mapping between FP and Schrödinger equation in terms of so-called Lévy–Schrödinger semigroup^31,32. Of course, direct numerical solutions are always possible.

The components (shown near curves) of oscillator strength (dipole matrix element Eq. (32)) for transition between the states 1s and 1p as functions of Lévy index \(\mu \). Both functions are strongly decreasing as \(\mu \rightarrow 0\).

Full size image

In this context, we also mention one more interesting physical problem regarding the onset of chaos³³ in the systems (like semiconductors^1,2 and trapped cold atomic gases³⁴) with the joint effect of the Lorentz force, Zeeman splitting, and spin-orbit coupling^35,36. This suggests one more generalization of our 3D fractional quantum-mechanical oscillator problem. Namely, the spin-orbit interaction term can be added to the fractional Schrödinger equation(4). In this case, the solution will be more sophisticated as the wave function will be spinor now (see, e.g.³⁷) although the problem can be solved in the momentum space similar to the present case. This problem turns out to be very important for amorphous semiconductors with atomic vibrations^1,2,38, where chaos can adversely influence possible optoelectronic, spintronic, and/or photovoltaic devices functionality.

Many physical problems can be reduced to that of a quantum harmonic oscillator. One of them is related to the theoretical studies of quark-antiquark bound states (so-called quarkonium, see, e.g. Ref.³⁹), which had been a hot topic several decades ago. In the context of 1D fractional quantum oscillator, the problem had been introduced by Laskin¹⁷, where the 1D potentials of the form \(|x|^\gamma \), \(1<\gamma \le 2\) had been considered instead of quadratic one. At the same time, as the consideration was phenomenological by its nature, the linear potential had also been applied for that problem³⁹. As it is well-known^18,19, such potential in the 1D Schrödinger equation admits exact solution in terms of Airy functions¹⁴. In our context, such problem arises in Eq. (11) for \(l=0\) and \(\mu =1\). This problem had been considered by us earlier⁴⁰ in the context of Lévy flights confinement by linear potential. The energies and wave functions of s-states (i.e. those for \(l=0\)) are in a very good agreement with those obtained from exact solutions in terms of Airy functions. A quantitative comparison is out of the frames of the present work. Note also, that 1D and 2D problems for the fractional Schrödinger equation in a linear potential had been considered in the paper⁴¹. For 1D case this paper complements both our studies of 1D fractional oscillator²⁰ for the entire \(\mu \) domain and those in the paper⁴⁰ for \(\mu =1\), where the exact solution in terms of Airy functions had been presented. The consideration of Cauchy distribution time evolution, made in⁴¹ is also complementary to our study⁴².

One more promising approach to non-relativistic quantum problems in central force potentials like \(r^q\), is to suggest an approximation to the spectrum of the problem for arbitrary q. This has been done in the paper⁴³ for the case of ordinary (i.e. non-fractional) laplacians in the Schrödinger equation. Moreover, the author⁴³ was instrumental to generalize the approximation both to the case of logarithmic potential, corresponding to \((dr^q/dq)_{q=0}=\ln r\) and negative q like hydrogenic problems with \(q=-1\). It would be interesting in future to generalize the approach⁴³ to the case of fractional laplacians in corresponding Schrödinger equations.

Several one-dimensional problems (like fractional quantum harmonic oscillator and infinite quantum well as a limiting case of the potentials \(|x|^p\) as \(p \rightarrow \infty \)) had been considered in Ref.⁴⁴. There, the corresponding time-dependent Schrödinger equation had been solved by the matrix method. Different nonlocal operators, based on fractional Riesz derivatives had been discussed at length. Numerical and variational solutions had been presented. The variational solution of 1D fractional quantum oscillator problem had been reported in terms of Hermite and Laguerre polynomial type functions. This approach is complementary to ours in the sense that here we considered the “Gaussian” trial functions for comparison. The approach of the paper⁴⁴ is also complementary to our 1D case, considered in Ref.²⁰. The consideration of the infinite potential well problem is complementary to our work⁴⁵, where the same problem have been solved by the explicit matrix method. Note that the matrix method had been applied in our numerical solutions of 3D (present consideration) and 1D fractional oscillators problems.

In summary, we have studied the spectral problem for a 3D fractional quantum harmonic oscillator with arbitrary Lévy index \(0<\mu \le 2\). Our main supposition here is that Laskin’s construction of path integrals with Lévy measure is equivalent to “extraction” of probability density function from stochastic fractional Langevin equation and, in turn, to the assumptions made in seminal Anderson paper²⁶. This permits us to assert that fractional Schrödinger equation accounts for disorder phenomenologically with Lévy index \(\mu \) being the measure of the degree of disorder. The presence of potential in the system “tames” the initial Lévy distribution, making it decay faster than that in a corresponding free problem. In other words, here once more we have an interplay between the width of disorder distribution and system potential, which makes the square of the modulus of the corresponding wave function decay faster in space. This makes the problem of a fractional 3D oscillator resemblant to Anderson localization²⁶ in disordered systems.

1. Snaith, H. J. Perovskites: The emergence of a new era for low-cost, high-efficiency solar cells. J. Phys. Chem. Lett. 4, 3623 (2013).

   Article  CAS  Google Scholar 

2. Stranks, S. D. & Snaith, H. J. Metal-halideperovskites for photovoltaic and light-emitting devices. Nat. Nanotechnol. 10, 391 (2015).

   Article  ADS  CAS  Google Scholar 

3. Zhu, H. et al. Lead halide perovskite nanowire lasers with low lasing thresholds and high quality factors. Nat. Mater. 14, 636 (2015).

   Article  ADS  CAS  Google Scholar 

4. Xing, G. et al. Low-temperature solution-processed wavelength-tunable perovskites for lasing. Nat. Mater. 13, 476 (2014).

   Article  ADS  CAS  Google Scholar 

5. Dong, S., Liu, J.-M., Cheong, S.-W. & Ren, Z. Multiferroic materials and magnetoelectric physics: Symmetry, entanglement, excitation, and topology. Adv. Phys. 64, 519 (2015).

   Article  ADS  CAS  Google Scholar 

6. Samko, S. G., Kilbas, A. A. & Maritchev, O. I. Fractional Integrals and Derivatives (Gordon and Breach, 2003).

   Google Scholar 

7. Podlubny, I. Fractional Differential Equations (Academic Press, 1999).

   MATH  Google Scholar 

8. Metzler, R. & Klafter, J. The random walk’s guide to anomalous diffusion a fractional dynamics approach. Phys. Rep. 339, 1 (2000).

   Article  ADS  MathSciNet  CAS  Google Scholar 

9. Metzler, R. & Klafter, J. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A Math. Gen. 37, R161 (2004).

   Article  ADS  MathSciNet  Google Scholar 

10. Lévy, P. Théorie de l’Addition des Variables aléatoires (Gauthier-Villars, 1954).

    MATH  Google Scholar 

11. Lévy Flights and related topics in physics (eds. Shlesinger, M. F., Zaslavsky, G. M. & Frisch, U. ). in Lecture Notes in Physics. (Springer, 1995).

12. Hughes, B. D. Random Walks and Random Environments Vol. 1 (Clarendon Press, 1995).

    MATH  Google Scholar 

13. Samorodnitsky, G. & Taqqu, M. S. Stable Non-Gaussian Random Processes (Chapman and Hall, 1994).

    MATH  Google Scholar 

14. Handbook of Mathematical Functions (eds. Abramowitz, M. & Stegun, I.) (Dover, 1972).

15. Jespersen, S., Metzler, R. & Fogedby, H. C. Lévy flights in external force fields: Langevin and fractional Fokker-Planck equations and their solutions. Phys. Rev. E 59, 2736 (1999).

    Article  ADS  CAS  Google Scholar 

16. Laskin, N. Fractional quantum mechanics. Phys. Rev. E 62, 3135 (2000).

    Article  ADS  CAS  Google Scholar 

17. Laskin, N. Fractional Quantum Mechanics (World Scientific, 2018).

    Book  Google Scholar 

18. Landau, L. D. & Lifshits, E. M. Quantum Mechanics. Nonrelativistic Theory (Pergamon Press, 1995).

    Google Scholar 

19. Flügge, S. Practical Quantum Mechanics (Springer, 1999).

    MATH  Google Scholar 

20. Kirichenko, E. V. & Stephanovich, V. A. Confinement of Lévy flights in a parabolic potential and fractional quantum oscillator. Phys. Rev. E 98, 052127 (2018).

    Article  ADS  CAS  Google Scholar 

21. Polyanin, A. D. & Zaitsev, V. F. Handbook of Exact Solutions of Ordinary Differential Equations (CRC Press, 1995).

    MATH  Google Scholar 

22. Kirkwood, J. R. An Introduction to Analysis (PWS Publishing, 1995).

    MATH  Google Scholar 

23. Kirichenko, E. V. & Stephanovich, V. A. The influence of Coulomb interaction screening on the excitons in disordered two-dimensional insulators. Sci. Rep. 11, 11956. https://doi.org/10.1038/s41598-021-91414-w (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

24. Kirichenko, E. V. & Stephanovich, V. A. The influence of disorder on the exciton spectra in two-dimensional structures. Phys. Chem. Chem. Phys. 21, 21847 (2019).

    Article  CAS  Google Scholar 

25. Stephanovich, V. A. & Olchawa, W. Lévy distributions and disorder in excitonic spectra. Phys. Chem. Chem. Phys. 22, 24462 (2020).

    Article  CAS  Google Scholar 

26. Anderson, P. W. Absence of diffusion in certain random lattices. Phys. Rev. 109, 1492 (1958).

    Article  ADS  CAS  Google Scholar 

27. Ashkroft, N. W. & Mermin, N. D. Solid State Physics (Harcourt, 1976).

    Google Scholar 

28. Laguta, V. V., Stephanovich, V. A., Raevski, I. P., Raevskaya, S. I., Titov, V.V., Smotrakov, V.G., & Eremkin, V.V. Magnetoelectric effect in antiferromagnetic multiferroic Pb\(({{\rm Fe}}_{1/2}{{\rm Nb}}_{1/2}){{\rm O}}_3\) and its solid solutions with \({{\rm PbTiO}}_3\). Phys. Rev. B 95, 014207 (2017).

29. Stephanovich, V. A. & Laguta, V. V. Transversal spin freezing and re-entrant spin glass phases in chemically disordered Fe-containing perovskite multiferroics. Phys. Chem. Chem. Phys. 18, 7229 (2016).

    Article  CAS  Google Scholar 

30. Zagorodniy, Yu. O., Kuzian, R. O., Kondakova, I. V., Maryško, M., Chlan, V., Štěpánková, H., Olekhnovich, N. M., Pushkarev, A. V., Radyush, Yu. V., Raevski, I. P., Zalar, B., Laguta, V. V., & Stephanovich, V. A. Chemical disorder and \(^{207}\)Pb hyperfine fields in the magnetoelectric multiferroic Pb\(({\text{Fe}}_{1/2}{\text{ Sb }}_{1/2}){\text{ O }}_3\) and its solid solution with Pb\(({\text{ Fe }}_{1/2}{\text{ Nb }}_{1/2}){\text{ O }}_3\). Phys. Rev. Mater. 2, 014401 (2018).

31. Risken, H. The Fokker-Planck Equation (Springer, 1996).

    Book  Google Scholar 

32. Pavliotis, G. A. Stochastic Processes and Applications (Springer, 2014).

    Book  Google Scholar 

33. Reichl, L. E. The Transition to Chaos. Conservative Classical Systems and Quantum Manifestations Vol. 2 (Springer, 2004).

    MATH  Google Scholar 

34. Larson, J., Anderson, B. M. & Altland, A. Chaos-driven dynamics in spin-orbit-coupled atomic gases. Phys. Rev. A 87, 013624 (2013).

    Article  ADS  Google Scholar 

35. Bychkov, Yu. A. & Rashba, E. I. Properties of a 2D electron gas with lifted spectral degeneracy. JETP Lett. 39, 78 (1984).

    ADS  Google Scholar 

36. Kirichenko, E. V., Stephanovich, V. A. & Sherman, E. Y. Chaotic cyclotron and Hall trajectories due to spin-orbit coupling. Ann. Phys. (Berlin) 2000012, 1–8. https://doi.org/10.1002/andp.202000012 (2020).

    Article  CAS  Google Scholar 

37. Stephanovich, V.A., Sherman, E. Ya., Zinner, N.T., &Marchukov, O.V. Energy-level repulsion by spin-orbit coupling in two-dimensional Rydberg excitons. Phys. Rev. B 97, 205407 (2018).

38. Fu, M. et al. Neutral and charged exciton fine structure in single lead halide perovskite nanocrystals revealed by magneto-optical spectroscopy. Nano Lett. 17, 2895 (2017).

    Article  ADS  CAS  Google Scholar 

39. Berkelman, K. Heavy quark physics. Rep. Prog. Phys. 49, 1 (1986).

    Article  ADS  CAS  Google Scholar 

40. Garbaczewski, P. & Stephanovich, V. Lévy flights in inhomogeneous environments. Physica A 389, 4419 (2010).

    Article  ADS  MathSciNet  CAS  Google Scholar 

41. Liemert, A. & Kienle, A. Fractional Schrödinger equation in the presence of the linear potential. Mathematics 4, 31. https://doi.org/10.3390/math4020031 (2016).

    Article  MATH  Google Scholar 

42. Garbaczewski, P. & Stephanovich, V. Lévy flights and nonlocal quantum dynamics. J. Math. Phys. 54, 072103 (2013).

    Article  ADS  MathSciNet  Google Scholar 

43. Hall, R. L. A simple interpolation formula for the spectra of power-law and log potentials. J. Phys. G Nucl. Part. Phys. 26, 981 (2000).

    Article  ADS  CAS  Google Scholar 

44. Herrmann, R. Solutions of the fractional Schrödinger equation via diagonalization—A plea for the harmonic oscillator basis. Part 1: The one dimensional case. arXiv:1805.03019 (2018).

45. Kirichenko, E. V., Garbaczewski, P., Stephanovich, V. & Żaba, M. Lévy flights in an infinite potential well as a hypersingular Fredholm problem. Phys. Rev. E 93, 052110 (2016).

    Article  ADS  Google Scholar 

Download references

Authors and Affiliations

1. Institute of Physics, University of Opole, ul. Oleska 48, 45-052, Opole, Poland

   V. A. Stephanovich & E. V. Kirichenko

2. Department of Physics and Medical Engineering, Rzeszów University of Technology, al. Powstańców Warszawy 6, 35-959, Rzeszów, Poland

   V. K. Dugaev

3. Department of Mathematics, Universidad de Las Palmas de Gran Canaria, Campus de Tafira Baja, C.P. 35017, Las Palmas, Spain

   Jackie Harjani Sauco & Belén López Brito

Contributions

V.S., V.D. and B.L.B. conceived the idea of the article. E.K., V.S., V.D. and J.H.S. wrote the main manuscript text. V.S., V.D. and B.L.B. performed analytical calculations and E.K., J.H.S. numerical ones. E.K. prepared Fig. 1, E.K. and J.H.S. prepared Fig. 2; V.D. and B.L.B. prepared Fig. 3. All authors reviewed the manuscript.

Corresponding author

Correspondence to V. A. Stephanovich.

Competing interests

The authors declare no competing interests.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions