Anomalous Relaxation and Three-Level System: A Fractional Schrödinger Equation Approach

: We investigate a three-level system in the context of the fractional Schrödinger equation by considering fractional differential operators in time and space, which promote anomalous relaxations and spreading of the wave packet. We ﬁrst consider the three-level system omitting the kinetic term, i.e., taking into account only the transition among the levels, to analyze the effect of the fractional time derivative. Afterward, we incorporate a kinetic term and the fractional derivative in space to analyze simultaneous wave packet transition and spreading among the levels. For these cases, we obtain analytical and numerical solutions. Our results show a wide variety of behaviors connected to the fractional operators, such as the non-conservation of probability and the anomalous spread of the wave packet.


Introduction
Many phenomena in nature, such as in complex viscoelastic media [1,2], electrical spectroscopy impedance [3][4][5], wave propagation in porous media [6,7], microflows of viscoelastic fluids [8], and gas transport in heterogeneous media [9][10][11] are associated with different types of scenarios, such as non-Markovian processes, which can be connected to memory effects, long-range correlations, and long-range interactions, among others. In this context, fractional calculus brings new possibilities in which the standard approach is unsuitable for describing these effects in connection with the experimental data. With a simple extension of the differential operator to the non-integer index, significant changes are incorporated into the standard approach, such as memory effects, long-range correlations, and other effects in complex systems [12,13]. Fractional extensions have been shown to be more suitable for capturing the experimental behavior when compared with integer-order models [14,15]. For instance, extending the dengue epidemic model using the Caputo derivative is more accurate than the standard model [16]. In the photothermal phenomena, the fractional operators are considered an extension of Fourier's law to study the heat propagation in non-homogeneous solids, which present different regimes of diffusivity [17]. These different regimes of diffusivity, such as subdiffusion and superdiffusion, can be modeled with different values of fractional order. Another application of fractional calculus is in the study of cardiac rhythm [18]. The normal rhythms of the heart are associated with a fractional order close to one; however, signals generated by hypertensive patients are modeled by a different fractional order. The hypertension model with fractional extension also was studied in Ref. [19]. Other applications of fractional calculus are found in Refs. [20][21][22][23][24][25][26] and the references therein. These fractional operators, and the processes connected to them, have been widely analyzed in several fields, particularly in quantum mechanics, as a fundamental physics theory that describes the physical properties of nature in lighting on a subatomic and atomic level. The pioneer works of N. Laskin [27][28][29] (based on extensions of the path integral formulation [30]), lead us to a fractional Schrödinger equation and have been followed by other extensions incorporating fractional differential operators in time and space [31], non-local terms [32,33], constraints among the different spatial coordinates (comb-model) [34][35][36][37][38], and others [39,40]. Space and time fractional derivatives in the Schrödinger equation were considered to obtain soliton-type solutions in the optical field [41], while numerical solutions for the fractional Schrödinger equations were obtained in Ref [42]. Other applications can be found in Refs. [43][44][45][46][47]. It is important to stress that, in these scenarios, anomalous relaxation processes [48] are present due to the extension of the differential operators of integer order to fractional ones. It is also worth mentioning that, from the analytical and numerical point of view, it is a challenge to obtain solutions when fractional derivatives in time and space are considered.
The main goal is to investigate the implications of considering time-dependent potentials in the fractional Schrödinger equation [49] given by where the fractional differential operator is the Caputo fractional time derivative, written as where k − 1 < α < k, ψ (k) ( r, τ) = ∂ψ (k) ( r, τ)/∂τ k , and k = α [31]. The fractional time derivative is defined as a convolution integral between a kernel (power law) and the wave function. We apply analytical and numerical methods to analyze Equation (1). As mentioned in Ref. [49], we can also extend the Schrödinger equation as follows: ih α ∂ α ∂t α ψ( r, t) = H(t)ψ( r, t).
Equations (1) and (3) are possible extensions and have been considered in many situations, such as optical solitons [50], resonant mode conversions [51], and two-level systems [52] as an extension of Rabi's model. However, note that when a Wick rotation is performed, the imaginary unit is increased to the same power as the time coordinates considered for Equation (1). Another point between the two equations involves the temporal behavior of the solution, which for the first case, is more suitable than for the second one, which decreases or grows with time instead of demonstrating a sinusoidal behavior. For these reasons, we consider Equation (1) in our developments.
Hence, we consider a three-level system in Equation (1) perturbed by an external electromagnetic field. Firstly, we analyze the system described by the time-dependent Hamiltonian in the form where the terms γ ij e ±iω ij t stand for the external perturbation, here understood as an external oscillatory field, with a frequency equal to ω ij . The terms γ ij are associated with the coupling between i and j, and E k means the energy eigenvalues of |k with k = 1, 2, 3 (E 2,3 > E 1 ). After analyzing the dynamics of the Schrödinger equations under the Hamiltonian given by Equation (4), we incorporate the kinetic terms with fractional powers yielding Note that the kinetic terms in Equation (5) can be related to a spatial fractional derivative, i.e., , refer to the Fourier transforms, respectively. This definition corresponds to the Riesz derivative [53,54].
We consider a three-level system because the atom shows a wide variety of effects in its interaction with the electromagnetic field, such as laser cooling, lasing without inversion, loss-free pulse propagation, and population transfer, compared to the two-level atom. Jing and Yu [55] presented an accurate quantum trajectory approach for a dissipative three-level model and developed a convolutionless stochastic Schrödinger equation called the time-local quantum state diffusion (QSD) equation without Markov approximation. They showed that this exact time-local QSD equation is suitable for exploring quantum dynamics for a higher dimensional quantum system with a non-Markovian environment. Castanos et al. [56] showed that the atoms with Λ and Ξ configurations present double points. This means that there are two independent quantum states with the same energy, which depends on the values of the dipole strengths of the interaction between the atoms and the radiation field. Robustness against the noise source is a big challenge for realizing quantum devices. In this way, based on generic three-level systems, Petiziol et al. [57] proposed a protocol according to the widespread stimulated Raman adiabatic passage technique. Their protocol realizes adiabatic acceleration by applying additional control fields to the optical excitation. They illustrated that the proposed method is efficient over various control parameters and brings the timescales close to the quantum speed limit when an environmental disturbance exists. Other applications of three-level systems can be found in Refs. [58][59][60] and the references therein.
Due to the additional degree of freedom, three-level systems are richer dynamically when compared with two-level ones. However, little attention has been devoted to them, mainly in the context of fractional quantum mechanics. In this context, we dedicate this research to studying a three-level system under the influence of an external perturbation. Our contribution to the fractional quantum mechanics field is to propose analytical and numerical solutions for this system. Aiming to understand the effects of fractional operators, we first analyze the Hamiltonian in the context of standard quantum mechanics, i.e., with integer-order operators. After, we extend the differential operators to a non-integer order by incorporating the fractional differential operators in that time and space. Our results suggest that the probability density follows an anomalous relaxation process due to the fractional operators. These developments are performed in Sections 2 and 3 by considering the analytical and numerical calculations for these Hamiltonians and analyzing the spreading behavior of the wave packet in different conditions. Finally, we present our discussions and conclusions in Section 4.

Schrödinger Equation without Kinetic Terms
The standard Schrödinger equation is defined by where i is the imaginary unit, ψ( r, t) is the wave function,h is the Planck's constant, and H is the Hamiltonian [61]. The following Hamiltonian can give a three-level system Considering Equation (6) with the Hamiltonian given by Equation (7) and, for simplicity,h = 1, we obtain the following equations The system described by Equation (10) corresponds to a coupled three-level system, respectively. A particular case was considered in Ref. [52] for a two-level system. This system of linear differential equations can be analyzed from analytical and numerical points of view. From the analytical point of view, we can use the Laplace transform to obtain the solutions for a given initial condition. We apply this approach combined with the Fourier transform to obtain the solution when the kinetic term is present in the Hamiltonian, i.e., when the Hamiltonian is given by Equation (26). Applying the Laplace transform L{ψ i (t); s} = ψ i (s) and L −1 {ψ i (s); t} = ψ i (t), we achieve the system of algebraic equations written as Solving Equations (12) and (13) for the ψ k , we find that with 2 23 and ψ 1 (0) = 1, while the solutions for ψ 2 (s) and ψ 3 (s) read as and A particular case from the previous one arises when ω 12 = ω 23 = 0 with the corresponding solutions ψ k (s): It is worth mentioning that the inverse Laplace transform can be found using the method of partial fractions, which implies finding the roots of a third-degree polynomial. For example, the case E 1 = 1, E 2 = 2, E 3 = 3 with γ 12 = γ 23 = 1 results in the following equations By applying the inverse of the Laplace transform to the functions ψ k (s), we obtain the corresponding ψ k (t) given by and We can use the numerical approach presented in Ref. [62] to find the solution to the previous set of equations. Numerical solutions for the probability densities |ψ 1 | 2 (red line), |ψ 2 | 2 (green line), and |ψ 3 | 2 (blue line) related to Equations (8)-(10) are presented in Figure 1.
. The choice of parameters is arbitrary. We chose these values due to the fact that some behaviors are highlighted with them. The probability of finding the system in one of the states obeys an oscillatory behavior. When no frequency is considered (panel (a)), the system stays more intense in the state ψ 1 . However, when the oscillatory dynamics of the field are considered (ω ij = 1), the system transit between the first (ψ 1 ) and third states (ψ 3 ) using the second state (ψ 2 ) as an edge of transition. This transition shows that the ψ 2 is the coupling between ψ 1 and ψ 3 . When no frequency is considered (panel (a)), the system stays more intense in the state ψ 1 . However, when the oscillatory dynamics of the field are considered (ω ij = 1), the system transitions between the first (ψ 1 ) and third state (ψ 3 ) using the second state (ψ 2 ) as an edge of transition. This transition shows that the ψ 2 is the coupling between ψ 1 and ψ 3 . These solutions show that if particles populated the first state initially, the transition for the third state depends on the second one. Moreover, the transmission of the particles depends on the parameters of the coupling and external field properties, such as the frequency. We Now, we extend the previous scenario by incorporating fractional time derivatives. For this case, we have that Interestingly, the Caputo fractional time derivative produces an anomalous relaxation for the solutions, which are usually asymptotically characterized by power laws. Similar to the previous case, we use a numerical calculation to obtain the behavior of the solutions of this set of coupled equations for an initial condition. In particular, for determining the numerical solutions for Equations (23)-(25), we use the Predictor-Evaluate-Corrector-Predictor method (PECE) [63] (see Appendix A for more details). Considering = 0, and α = 0.9, the solutions for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 are displayed in Figure 2 by the red, green, and blue lines, respectively. Panel (a) exhibits the case with ω ij = 0, i.e., a constant potential. This situation shows how drastically the dynamics change front fractional operators. Now occurs a competition between the first and second states in terms of their probability amplitudes, while the third state is less populated. On the other hand, if the potential oscillates, i.e., ω ij = 1, the second state works just as an edge between the first and third states, which are populated with the same amplitude, as seen in Figure 2b. In both cases, the probability amplitude goes proportionally with 1/α 2 . This result is expected when considering i α , as found in Ref. [64] for a particular case. The resulting probability increases can be interpreted as the creation of particles. We

Schrödinger Equation with Kinetic Dependences
The kinetic terms can be considered in Equation (4) by the following transformation For simplicity, we consider m = 1 andh = 1 and the one-dimensional case. Taking into account the kinetic terms, the dynamics of the three-level system are described by The analytical solution for this set of equations can be obtained using the Laplace and Fourier (F {ψ i (x, t); k} = ψ i (k, t) and F −1 { ψ i (k, t); x} = ψ i (x, t)) transforms. Applying the integral transforms to the previous set of equations assuming the initial conditions ψ 1 (x, 0) = ϕ(x), ψ 2 (x, 0) = 0, and ψ 3 (x, 0) = 0, we obtain the corresponding set of algebraic equations, written in the form: From Equation (32), we relate the solutions ψ 3 (p, s) and ψ 2 (p, s) by means of The former result combined with Equation (31) fixes the function ψ 2 (p, s) in terms of ψ 1 (p, s + iω 12 ) according to which, applied to Equation (30), leads to By applying the inverse of the Laplace transform to Equation (35), we obtain the wave function ψ 1 (p, t) in the form where Similarly, the functions ψ 2 (p, t) and ψ 3 (p, t) read as with Λ 2 (t) and Λ 3 (t) determined from and The inverse of the Fourier transform of the previous equations is given by for i = 1, 2, and 3 with The discretization of Equations (27)- (29) can be obtained by the finite difference scheme [65]. For this, we consider a grid defined by [0, L] × [0, T] with boundary conditions equal to ψ 1,2,3 (0, t) = ψ 1,2,3 (L, t) = 0. The time is discretized according to t m = mh t , where m = 1, 2, ..., N t with a time step equal to h t = T/N t ; each space coordinate is given by x n = nh x where n = 1, 2, ..., N x with a space step equal to h x = L/N x (see Appendix B for more details). To avoid numerical boundary problems, the origin of our space coordinate is located in L/2. With these considerations, Equations (27)-(29) are given by The numerical solutions for Equations (43)- (45) are displayed in Figure 3 in panels (a)-(f). Panels (a)-(c) show the probability amplitude (color scale) distribution in time and space for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 . As time advances, the Gaussian package that is centered in |ψ 1 | 2 spreads along the second (ψ 2 ) and third states (ψ 3 ). The packet scattering is more evident when we look for the density plot in Figures 3d-f, for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 2 | 3 . These results show that the probability of finding particles in a given state decreases in amplitude as time advances and becomes wider in the space domain. Furthermore, the particles transition from the first state to the second one, and, just after populating the second one, they go into the third state. This oscillatory transition reflects the effects of the potential and external perturbation. For these results, we consider h t = 0.0001, h x = 0.2, ω 12 = ω 23 = π/2, γ 12 = γ 23 = 2, and the initial condition Figure 3. Probability amplitude (color scale) spread in panels (a-c) for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 . Density plot of the probability distributions in panels (d-f) for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 . We consider h t = 0.0001, h x = 0.2, h t /2h 2 x = 0.01, ω 12 = ω 23 = π/2, γ 12 = γ 23 = 2, and the initial condition The profiles of |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 are exhibited in Figure 4, respectively, by panels (a)-(c). The red line is for t = 0.1, the green line is for t = 1.8 in panel (a) and t = 1.6 in panels (b) and (c), and the blue line is for t = 3.8 in panel (a) and t = 10 in panels (b) and (c). From these profiles, we note the distribution of the Gaussian package among the three states. The Gaussian shape spreads along space as time evolves. (c) Figure 4. Profiles for the Gaussian package distribution for the probabilities |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 in panels (a-c), respectively. The red line is for t = 0.1, the green line is for t = 1.8 in (a) and t = 1.6 in (b,c), and the blue line is for t = 3.8 in (a) and t = 10 in (b,c). We consider h t = 0.0001, h x = 0.2, h t /2h 2 x = 0.01, ω 12 = ω 23 = π/2, γ 12 = γ 23 = 2, and the initial condition One way to measure a probability package's spread is by computing the mean square displacement ( (∆x) 2 ), which is displayed in Figure 5. The red line is associated with |ψ 1 | 2 , while the black and blue lines are associated with |ψ 2 | 2 and |ψ 3 | 2 , respectively. In dealing with the standard formulation, the displacements evolve proportionally linearly with t 2 , which characterizes the standard spreading of the wave packet. Due to the numerical imprecision and the potential effects, our respective slopes for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 are 1.95, 2.01, and 1.83. The oscillations in the displacements reflect the oscillatory character of the potential, which causes a non-linear transmission among the states. It is worth mentioning that the position in relation to the (∆x) 2 -axis is translated to consider all the curves in the same figure without superposition. This shift does not alter the slope, which is an important parameter for our study.
Let us extend the previous result by incorporating a fractional differential in space. For this case, the previous set of equations can be written as follows: The solution for the previous set of equations is given by for i = 1, 2, and 3 with which resembles the form of the Lévy distribution found in anomalous diffusion processes. In Equation (50) In Equation (51), the function χ(s) reads which involves Mellin-Barnes integrals [31] (for more details see also Ref. [66]). The asymptotic behavior of Equation (50) in the limit of |x| → ∞ is given by G µ (x, t) ∼ i t/ 2|x| 1+µ , which is different from the usual one characterized by the Gaussian behavior. Note that this result for the asymptotic limit can be obtained by using the approach employed in Ref. [67]. It is essentially an integration over the Mellin-Branes integral poles, which represents Equation (50). This feature is directly connected to the presence of spatial fractional differential operators in Equations (46)- (48). For numerical investigation, we consider the same method as in the standard case, i.e., the finite differences. Considering this scheme, the discretization of Equations (46)-(48), yields The inclusion of fractional space derivatives (e.g., µ = 1.5) in the Schrödinger threelevel systems implies a delay in the probability spread, as shown by Figure 6a-c for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 , respectively. Furthermore, the scattering gains more evidence of the potential, which is seen through a sinuous scattering, as displayed by panels (d)-(f), in Figure 6, for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 , respectively. Figure 6. Effects of fractional space derivative with µ = 1.5 in the probability amplitude (color scale) spread (panels (a-c), for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 ) and in the density plot of the probability distributions in panels (d-f) for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 . We consider h t = 0.0001, h x = 0.2, ω 12 = ω 23 = π/2, γ 12 = γ 23 = 2, and the initial condition Profiles for |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 are displayed in Figure 7, respectively, in panels (a-c). We note from the profiles that the spread has a higher amplitude than the standard case and occurs in sinusoidal form due to the potential effects. The higher amplitude shows the package spread with a low velocity when µ = 2 was considered. . Profiles for the Gaussian package distribution for the probabilities |ψ 1 | 2 , |ψ 2 | 2 , and |ψ 3 | 2 in panels (a-c), respectively. We consider µ = 1.5. The red line is for t = 0.1, the green line is for t = 1.8 in (a) and t = 1.6 in (b,c), and the blue line is for t = 3.8 in (a) and t = 10 in (b,c). We consider h t = 0.0001, h x = 0.2, ω 12 = ω 23 = π/2, γ 12 = γ 23 = 2, and the initial condition The mean square displacement for |ψ 1 | 2 (red line), |ψ 2 | 2 (black line), and |ψ 3 | 2 (blue line) with µ = 1.5 is exhibited in Figure 8. The mean square displacement increases by t 1.32 , for |ψ 1 | 2 ; t 1.28 , for |ψ 2 | 2 ; and t 1.27 , for |ψ 3 | 2 . The sinusoidal form reflects the potential effects and transitions among the states.

Discussions and Conclusions
We have investigated a three-level system described by a fractional Schrödinger equation, with a time-dependent potential and an external electromagnetic field. We considered fractional extensions of the Schrödinger equation in terms of differential operators of noninteger order in space and time. We first analyzed the three-level system by considering the transition among the levels to investigate how the probability evolves, in time, according to the conditions considered. Considering that one state is initially populated, our results suggest that the other two states become populated over time. This transition occurs in an oscillatory way, depending on the parameters related to the external field, e.g., oscillation frequency. We verified that extending the time operator to a fractional one introduces an anomalous relaxation and oscillation among the levels. Afterward, we incorporated a kinetic term into the previous scenarios. This feature allows the spreading of the wave packet at each level while the transition between the levels occurs. For this case, we also extended the differential operators to fractional differential operators in space. Extending the space operators to fractional ones, we verify an anomalous relaxation process related to the probability transition among the states. The transition process can be amplified according to the parameter's range. For example, if we address different weight coupling, the higher weight state will be populated earlier, and this state can bind the particles. Another way to change the probability transition is by changing the frequency. The frequency of the external field is directly connected with the probability transition frequency.
In summary, we have considered the fractional operator defined by a power-law kernel, which is well known as the Caputo derivative for the fractional time derivative. One of the most relevant aspects obtained with this extension of the time derivative is the possibility of describing different effects not present in the standard approach, such as memory effects and long-range correlations. These aspects have implications for the relaxation processes of these systems, which are different from the usual ones. These relaxation processes emerge from the fact that we may connect the solutions with the Mittag-Lefler functions, which are asymptotically governed by power laws. In future works, we will consider other fractional time derivatives with different kernels, which imply different relaxation processes. For the spatial operator, we considered the Riesz fractional operator, which also produces an anomalous relaxation with distributions asymptotically governed by power laws. Finally, we hope the results presented here can be useful in discussing anomalous processes and fractional Schrödinger equations.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Fractional Ordinary Differential Equations
One numerical method to solve fractional ordinary differential equations based on Caputo's definition is a generalization of the Adams-Bashforth-Moulton method, which is known as the Predictor-Evaluate-Corrector-Predictor method (PECE) [63].
Given an initial problem, defined by with y(0) = y 0 as the initial condition. The problem is solved, according to Ref. [63] by the following discretization The coefficients are defined by where j = 1, 2, ..., n, and n is associated with the discrete time window, T, which is discrete in t n = nh with n = 0, 1, ..., N where T = Nh.

Appendix B. Fractional Partial Differential Equations
In this work, Equations (46)- (48) are discretized according to the finite difference method [65]. The extensions for the fractional partial differential equations (FPDE) of this method are found in Refs. [68,69]. The first step in building this method is to consider a grid defined by [0, L] × [0, T]. The space and time are discretized by: x n = nh x , t m = mh t where n = 1, 2, ..., N x and m = 1, 2, ..., N t . In this way, the space and time steps are, respectively: h x = L/N x and h t = T/N t . To avoid numerical problems, we translate the origin of our space to L/2, without loss of generality. To solve the FPDE, it is necessary to be informed of the boundary conditions and initial conditions for all states. For simplicity, we consider boundary conditions equal to ψ 1,2,3 (0, t) = ψ 1,2,3 (L, t) = 0 and the initial conditions ψ 1 (x, 0) = e where J (n, k; µ) = (n − k) 2−µ − (n − k − 1) 2−µ . These equations are solved numerically by the Algorithm A1, described below.
We performed our numerical simulation using the language C, and the graphs were plotted using Gnuplot 5.2.

Appendix C. List of Symbols
A list of symbols used in this work is shown in Table A1.