Fractional quantum oscillator and disorder in the vibrational spectra.

We study the role of disorder in the vibration spectra of molecules and atoms in solids. This disorder may be described phenomenologically by a fractional generalization of ordinary quantum-mechanical oscillator problem. To be specific, this is accomplished by the introduction of a so-called fractional Laplacian (Riesz fractional derivative) to the Scrödinger equation with three-dimensional (3D) quadratic potential. To solve the obtained 3D spectral problem, we pass to the momentum space, where the problem simplifies greatly as fractional Laplacian becomes simply [Formula: see text], k is a modulus of the momentum vector and [Formula: see text] is Lévy index, characterizing the degree of disorder. In this case, [Formula: see text] corresponds to the strongest disorder, while [Formula: see text] to the weakest so that the case [Formula: see text] corresponds to "ordinary" (i.e. that without fractional derivatives) 3D quantum harmonic oscillator. Combining analytical (variational) and numerical methods, we have shown that in the fractional (disordered) 3D oscillator problem, the famous orbital momentum degeneracy is lifted so that its energy starts to depend on orbital quantum number l. These features can have a strong impact on the physical properties of many solids, ranging from multiferroics to oxide heterostructures, which, in turn, are usable in modern microelectronic devices.

Introduction

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[5]. 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[8-13], which are Markovian random processes whose probability density functions (pdf) are Lévy stable laws, characterized by the Lévy index . 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. ([8]). In dimensionless units (diffusion coefficient and particle mass are set to unity) it readswhere is an n-dimensional fractional Laplacian (Riesz fractional derivative[6,7])which at yields the ordinary one[6,7]. Here is -function[14], is (also n- dimensional) gradient operator and is an external potential, in which the anomalous diffusion occurs[15]. 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 . 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 (, where 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 space. The only difference is that the potential in this new equation will be . This implies that to solve our problem in momentum space, we can use well developed numerical and variational techniques.

The spectral problem for fractional quantum harmonic oscillator

In the units, where (m and are the mass and frequency of oscillating particle respectively), the spectral problem for 3D fractional oscillator readswhere the operator is defined by the expression (2) and is the eigenfunction of a fractional quantum harmonic oscillator having the eigenenergy for any specific value. In general 3D case, as it is typical for the quantum-mechanical problems in central force potentials, the eigenfunction 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 , and (here we suppress the indices for brevity) in momentum space render as second derivatives , and , the potential term in the equation (4) renders simply as an ordinary Laplacian in momentum space. At the same time, the fractional Laplacian in momentum space becomes simply , where is a vector modulus. Latter definition immediately implies that . This means that in momentum space, the equation (4) renders aswhere is the system Hamiltonian in momentum space. One note is in place here. Namely, the equation (5) shows that at , the famous orbital 3D oscillator degeneracy (see Ref.[18] 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 , the time-reversal symmetry remains (it is lifted, for instance, in the external magnetic field[18]), we still have degeneracy with respect to the quantum number m. Namely, the spectrum is invariant as . 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 admit the separation of angular and radial (in space) variables in the formwhere are spherical harmonics[14,18]. Here are associate Legendre polynomials of the first kind and is the corresponding normalization coefficient, which is obtained from the following normalization condition[14,18]where .

It is well-known from “ordinary” quantum mechanics[18,19] that 3D quantum oscillator problem for admits the factorization of the wave function (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 , and . In this case, the initial energy . After problem being in such a way factorized, the eigenenergy can be derived in the well-known form[18,19]where with being the quantum numbers, corresponding to respectively. Such factorization in the “ordinary” case 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 , and numbers.

It had been shown in[19], that the possibility of such factorization at is based on additivity of both Laplacian operator (kinetic energy) and the potential . At the same time, in 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 space, where the potential is factorable for (ordinary case) only.

This signifies the difference between fractional and ordinary 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 Here we also suppress the indices 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 functionfor which the Eq. (9) assumes particularly simple formwhich 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 . We note here, that if for the states with , the function R(k) has a maximum at , the corresponding function 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.

The Eqs. (9), (11) represent Schrödinger equation with the effective 1D potential . The plots of this potential at different and l are shown in Fig.1. It is seen from Fig.1a that for the problem is identical to that for 1D fractional oscillator[20], while for (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 as k is now the modulus of a momentum vector. As it is seen from Fig.1b, the centrifugal term in the potential is dominant at , while at the main contribution comes from . This means that the whole spectrum of 3D fractional oscillator remains to be discrete similar to the “ordinary” case . Note that at small the potential begins to grow at very large k, while at the potential goes to 1 at . The same situation at gives constant value one for , see Fig. 1a. Our numerical analysis shows that in 3D case despite the contribution of centrifugal term, at , the whole oscillator spectrum collapses into one value, see below. This property is similar to that in 1D case[20]. Note that at and the potential at the potential has nonanalytical behavior with infinite first derivative. At this asymptotics is “killed” by the centrifugal term.

Figure 1

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

Variational solution. Classification of the states of 3D fractional oscillator

It is well-known that the variational principle of quantum mechanics[18] 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 . Such “simple” trial functions can be constructed on the base of the Gaussian one, inherent in an “ordinary” (i.e. that at ) 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 , given by the radial (in 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 . This yields the equation for large k in the form . The spatially decaying solution to this equation is proportional to , where and [21]. Here is so-called MacDonald function with the following large x asymptotics [14]. Substitution of latter asymptotics into the above expression generates following wave function asymptotics at all admissible ’sHere we once more suppressed indices nlm. In the case of (ordinary quantum oscillator) the asymptotics (Eq. 12) reproduces well known result for “ordinary” 3D quantum oscillator[18]. It is seen that for all the asymptotics (Eq. 12) decays sufficiently fast for corresponding integrals to be convergent. Namely, in worst case of , 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.[22]. Our analysis show that at small k the functions R(k) have the same asymptotics (and accordingly ) as in ordinary case, corresponding to [18].

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 instead of 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 there is no obvious way to introduce the principal quantum number n. We recollect that for ordinary 3D oscillator (see, e.g.[19]), where 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 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 to this state. The next (excited) state for the same l has one node, i.e. etc. This means that we may classify the states of our fractional oscillator by essentially two quantum numbers - and l. To be specific, for (s - state; here we use the common notations of the states with the given l, see, e.g.[18]) we have the whole sequence ; the same is true for . Introducing the principal quantum number , we obtain following explicit classification of the states of the fractional quantum oscillator: 1s (, hence , ), 1p (, ), 1d (, ),... 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 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 and l. Specifically, in terms of functions , we haveHere, are normalization constants, found from the “1D-like” conditionsAlso, are variational parameters, while the constants , and (here ) are found from the orthogonality conditionswhere is Cronecker delta. We note that the trial functions for different l (for instance and ) 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 (and the same l) should be orthogonalized. Namely, in the particular case of the above trial functions we have explicitlyThe orthogonalization procedures Eqs. (17), (18) along with normalization of the functions determine uniquely the coefficients through variational parameters .

We have explicitly Here is the -function[14]. The other normalization coefficients can be calculated, but their explicit expressions become progressively more cumbersome.

Orthogonality condition yieldsSubstitution of the expression (20) into the normalization coefficient Eq. (19b) generates following explicit expression for the function where is determined by the expression (20) andNow, after variational determination of the parameter (see below), the parameter (20) starts to depend on the single parameter , which can also be determined variationally. Having determined the parameters and , we then use the orthogonality condition to express both and through . After substitution of the above found and into the corresponding normalization coefficient, the function starts to depend on only. Such “chain rule” can be applied to determine uniquely all higher order functions, corresponding to (s-state), (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 . Namely, the arbitrary trial function is the product of the factors (small k asymptotics), (large k asymptotics) and a polynomial of degree 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 energieswhere is the Hamiltonian (5) and are the variational approximations of the eigenenergies . Normally . As the radial parts 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 formwhere . 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 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 into (24) gives the approximate value of the ground state energy for all , the same with, say, gives the energy of the excited state 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 .

Substitution of the trial function into the integral (Eq. 24) with subsequent minimization over yieldswhich is similar to the expression for 1D fractional oscillator[20]. Further substitution of this value to the expression (24) generates the approximate value of the ground state energy for arbitrary It is seen that gives correct value 3/2 of the ground state energy for , corresponding to the ordinary quantum oscillator with the spectrum in our units. Below we shall see, that for the case all the spectrum shrinks into a single value , 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 gives that , which is given by Ex. (25). The variational expression for the energy of the 1p state readsThe expression (27) shows that at 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 , see above. As for 1p state we have (no nodes) and so that and the energy . At the same time, at we have removable divergence , while the rest of the expression (27) gives , 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 . Within the suggested variational approach, the energy of any desired level can be evaluated analytically, although the calculations for higher excited states become very cumbersome, see above.

Numerical analysis

The easiest way to solve our problem numerically is to use the Eq. (11) for the function , which has only one divergency (in the centrifugal term ) at . The numerical solution of the Eq. (11) should be augmented by the boundary conditions . 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 10000, which makes the task quite computer intensive.

As we can reproduce numerically the spectrum of the problem for arbitrary , 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 functions (10), we havewhere “G” in lower index stands for “Gaussian”, A is normalization coefficient and is variational parameter. We note that at the function is equivalent to (see (13)) with .

The normalization of (28) with the help of condition (16) yieldsSubstitution of Eq. 28 with respect to Eq. 29 into the expression for energy (Eq. 24) with its subsequent minimization over generates following minimizing wave function widthFurther 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 is determined by Eq. (30). At , the energy also gives correct value 3/2 of the oscillator ground state energy. The same is valid for , where . Here we also have a removable divergency of the type 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[18]. It is seen that the agreement is a little worse for the 2s state, which has one node. For the nodeless states with , 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 and , the relative error should be maximal somewhere inside this interval. The inset in Fig. 2c shows that the maximal relative error occurs at 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.

Figure 2

(a) Reports the lowest state eigenenergies (indicated near the curves) as functions of Lévy index , 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 domain. (d) Portrays the comparison between variational (both types of the functions, shown in the legend) and numerical 1s wave functions for the same values, as in (b). Insets report the relative errors for and 0.4 as the wave function for 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 .

The selected numerical wave functions in -space are reported in Fig. 2b,d for different values of . 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 have different decay rates, which are dictated by the asymptotics (Eq. 12). For instance, the 2s function for at is not yet achieved its final part (compare to that for ), which occurs at . For smaller , the spatial extension of the wave functions in space increases drastically. The functions for in space are almost delocalized, which means that in coordinate space they are Dirac -function-like. This behavior is qualitatively similar to that of 1D fractional quantum oscillator[20] as well as for 2D[23,24] and 3D[25] 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 , which is close to the above value , where the relative error in the energy is maximal. As numerical and variational curves are identical for (the corresponding trial functions are exact), the relative error is reported for 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 and 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 as for 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 space. Such good variational approximation by the functions, having only one parameter may be because the energies are monotonous functions of . 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.

One more observation is in place here. Namely, the eigenenergies in fractional case are smaller than those for “ordinary” case . 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.[23] 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[20], i.e. the slowest decay have the functions for (“ordinary” quantum oscillator), while the fastest (like Dirac , see above) have the functions for ( in Fig. 3), corresponding physically to the strongest disorder. This may be one more signature of the famous Anderson localization[26] phenomenon.

Figure 3

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

Possible experimental applications. Outlook

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 ) for the transition between the energy levels with different orbital quantum numbers l and . 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 and . As here we are working in momentum space, the formula for the oscillator strength contains the matrix element 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 and respectively. For the transition between the states 1p and 1s, the explicit expressions for matrix elements assume the form whereSubstitution of the matrix elements (Eqs. 33a–33c) with respect to numerically calculated and wave functions into the expression (32), generates the Lévy index dependence of oscillator strength components (), which is reported in Fig. 4. It is seen that the oscillator strength goes to zero in the case of strongest disorder, corresponding to , see Ref.[23] for details. We observe (see above) that as , both and go to the same value (as, actually, the whole spectrum) 1/2. This means that for , 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 -functions. This means that the probability of optical transitions is much higher in relatively ordered regions, where . 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[5]. The non-Gaussian statistics due to disorder and frustration plays an important role in these substances[28-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.

Figure 4

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 . Both functions are strongly decreasing as .

In this context, we also mention one more interesting physical problem regarding the onset of chaos[33] in the systems (like semiconductors[1,2] and trapped cold atomic gases[34]) 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.[37]) 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.[39]), which had been a hot topic several decades ago. In the context of 1D fractional quantum oscillator, the problem had been introduced by Laskin[17], where the 1D potentials of the form , 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[39]. As it is well-known[18,19], such potential in the 1D Schrödinger equation admits exact solution in terms of Airy functions[14]. In our context, such problem arises in Eq. (11) for and . This problem had been considered by us earlier[40] in the context of Lévy flights confinement by linear potential. The energies and wave functions of s-states (i.e. those for ) 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[41]. For 1D case this paper complements both our studies of 1D fractional oscillator[20] for the entire domain and those in the paper[40] for , where the exact solution in terms of Airy functions had been presented. The consideration of Cauchy distribution time evolution, made in[41] is also complementary to our study[42].

One more promising approach to non-relativistic quantum problems in central force potentials like , is to suggest an approximation to the spectrum of the problem for arbitrary q. This has been done in the paper[43] for the case of ordinary (i.e. non-fractional) laplacians in the Schrödinger equation. Moreover, the author[43] was instrumental to generalize the approximation both to the case of logarithmic potential, corresponding to and negative q like hydrogenic problems with . It would be interesting in future to generalize the approach[43] 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 as ) had been considered in Ref.[44]. 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[44] is also complementary to our 1D case, considered in Ref.[20]. The consideration of the infinite potential well problem is complementary to our work[45], 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 . 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[26]. This permits us to assert that fractional Schrödinger equation accounts for disorder phenomenologically with Lévy index 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[26] in disordered systems.

Methods

The details of our theoretical methodology and those of working with fractional derivatives and fractional Laplacians, in particular, have been described in the sections “The spectral problem for fractional quantum harmonic oscillator” and “Variational solution. Classification of the states of 3D fractional oscillator”. The numerical solutions of boundary value problems for partial differential equation (11) have been conducted using the commercial Mathematica software package.