Anomalous diffusion of the heavy quarks through the fractional Langevin equation

The dynamics of heavy quarks within the hot QCD medium have been revisited, considering the influence of anomalous diffusion. This study has been conducted using the framework of the fractional Langevin equation involving the Caputo fractional derivative. We introduce a numerical scheme for the fractional Langevin equation and demonstrate that the mean square displacement of the particle exhibits anomalous diffusion, deviating from a linear relationship with time. Our analysis calculates various entities, such as mean squared momentum, momentum spread, and the nuclear suppression factor, $R_{AA}$. Notably, our findings indicate that superdiffusion strongly suppresses the $R_{AA}$ compared to normal diffusion in the hot QCD medium. The possible impacts on other parameters are also discussed.


I. INTRODUCTION
The ultra-relativistic heavy-ion collisions (HICs) at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) have predicted the existence of Quark-Gluon Plasma (QGP); a state of matter where quarks and gluons are free to move beyond the nucleonic volume [1][2][3][4][5].The QGP is short-lived, with an expected lifespan of a few fm/c (approximately 4-5 fm/c at RHIC and 10-12 fm/c at LHC [6,7]).Conversely, investigating the characteristics of the QGP by studying the dynamics of heavy quarks (charm and beauty) is a subject of significant interest.In this context, The heavy quarks (HQs) serve as prominent probes of QGP .Due to their large masses (M c /M b ∼ 1.3/4.5 GeV), HQs are generated at the early stages of HICs.Moreover, the thermalization of the HQs is delayed compared to the light partons of the bulk medium by a factor on the order of ∼ M/T .This delay renders the thermalization time of the HQs comparable to the lifetime of the QGP fireball.Since the HQs are not expected to achieve complete thermalization, they retain the memory of their interaction history.Hence, they serve as a novel probe for observing the complete evolution of the QGP medium and act as non-equilibrium entities within the equilibrated QGP.The standard Langevin equation (LE) is used to study Brownian motion under the assumption of normal diffusion, which means a linear increase in the mean-square displacement of the particle with time for a sufficiently large time [38].By using the Langevin equation [39][40][41][42][43] in the context of the HQs momentum and position evolution in the hot QCD medium and numerous studies available for the experimental observables related to the HQs, such as the nuclear suppression factor, R AA and elliptic flow, v 2 [26,[44][45][46].However, the conventional approach faces challenges in adequately explaining the experimental observations concerning the energy loss experienced * jaiprakashaggrawal2@gmail.com by HQs within the QGP medium [46].This prompts the question: Could the anomalous motion of HQs offer a more suitable explanation for the observed data?Addressing this question necessitates an investigation into whether the anomalous behaviour observed in the motion of non-relativistic Brownian particles can also occur in the motion of relativistic Brownian particles.By examining this possibility, we can determine whether anomalous motion could potentially provide insights into the discrepancies observed in the energy loss of HQs within the QGP medium.Several processes exhibit anomalous diffusion, a phenomenon where the mean squared displacement of the particle does not vary linearly with time as predicted by Einstein's law.Specifically, anomalous diffusion is characterized by a non-standard time dependence of the mean squared displacement (⟨x(t) 2 ⟩) written as follows [47][48][49][50][51][52][53][54][55][56][57][58][59][60], where t is the evolution time of the particle.The process describes normal diffusion, corresponding to the case where µ = 1.Subdiffusion occurs when µ < 1, whereas superdiffusion occurs when µ > 1 [49].Despite the crucial role played by the LE across several fields, it fails to accurately describe certain behaviours, such as anomalous diffusion (i.e., superdiffusion and subdiffusion).Therefore, fractional Langevin equations (FLE) have been proposed to study anomalous diffusion [61][62][63][64].Mainardi et al., [65][66][67] introduced an FLE in their groundbreaking research.The Langevin method has primarily been applied to study anomalous diffusion using Caputo fractional derivative [68][69][70].On the other hand, there are several distinct formulations for fractional derivatives, such as Riemann-Liouville derivative [71], Riesz derivative [72], Feller derivative [73], and others.In the subsequent, we will only use the Caputo fractional derivative in our analysis to study the anomalous diffusion of the particle in the medium.Another valuable approach to studying anomalous diffusion involves the investigation of the fractional diffusion equation, fractional Fokker-Planck equation [74], and generalized Chapman-Kolmogorov equation [75].Anomalous diffusion has been used in many fields, including molecular chemistry [76], biology [77], and anomalous diffusion, polymer transport theory [78].Also, when accounting for interactions between the Brownian particle and the constituent particles of the medium, the fluctuations are influenced by their prior states, a phenomenon termed as memory.The current movement is affected by previous movements via a memory kernel in the generalized Langevin equation.Such memory effects can result in anomalous diffusion for a particular memory time [79][80][81][82][83][84][85] and also plays a role in the QGP [86][87][88][89][90].
Recently, the transverse momentum broadening of a fast parton via superdiffusion in the QCD matter has been studied in Ref. [91].
With this motivation for the applications of superdiffusion on the HQs dynamics in hot QCD matter, we use FLE, a generalized form of the LE.Unlike the LE, the FLE replaces integer-order derivatives with fractionalorder derivatives, specifically of Caputo type [92].This may be the first attempt where we present an anomalous diffusion for the HQ dynamics in the QGP medium.In this article, we specifically study the effect of superdiffusion on the HQ dynamics in the QGP medium.The subdiffusion is not considered within the scope of the current analysis.
We anticipate a key outcome in our study, specifically that in the presence of superdiffusion, there will be an increase in the energy loss experienced by the HQs within the QGP medium.This conclusion is confirmed by studying various key parameters, including the mean squared displacement, mean squared momentum, momentum spread, dN/dp T , and the R AA of the HQs.Notably, the latter holds particular significance in the phenomenology of the HQs within the QGP.The strong suppression observed in R AA is anticipated to contribute significantly to developing a substantial elliptic flow (v 2 ) for the HQs.
The article is structured as follows: Section II introduces the formalism; we analytically solve the FLE of non-relativistic heavy particles using the Laplace technique, along with presenting a numerical scheme for solving the relativistic FLE.Section III is dedicated to the presentation of our results.Lastly, Section IV provides a comprehensive summary of our conclusions.

II. THE FRACTION LANGEVIN EQUATION FOR BROWNIAN PARTICLES
For illustrative purposes, we discuss the scenario of a heavy particle with mass M undergoing one-dimensional (1-D) motion within the nonrelativistic limit.We use the FLE to describe the dynamics of this particle.In this framework, the evolution of the particle's position and momentum is characterized by fractional derivatives of orders β and α [61,62], where the momentum of a particle at time t is denoted as p(t), its position as x(t).C D α 0+ denotes Caputo fractional derivative, α and β are the fractional parameter, with n − 1 < α ≤ n, and C D β 0+ with n − 1 < β ≤ n, n ∈ N (n is a natural number).Brownian particle encounters two distinct forces: the dissipative force, characterized by the drag coefficient, denoted by γ, and the stochastic force, denoted as ξ(t).The latter governs the random noise, commonly referred to as white Gaussian noise.White noise gives rise to a fluctuating field without memory, characterized by instantaneous decay in correlations of the white noise, often referred to as a δ correlation.The random force satisfies certain properties, such as: D is the diffusion coefficient of the heavy particle in a medium of temperature, T .The drag coefficient (γ) is related to the diffusion coefficient through the Fluctuation-Dissipation Theorem (FDT) as follows: The Caputo fractional derivative [92] is defined as, where u (n) denotes the n th derivative of u.The fractional derivative which corresponds to the superdiffusion, i.e., when 1 < ν ≤ 2, is given by where Γ(•) denotes the gamma function.In the following subsection, we solve analytically the FLE of nonrelativistic Brownian particles.

A. Analytical solutions
Analytically, we solve the FLE governing the motion of non-relativistic Brownian particles.To achieve this, we use the Laplace technique, providing analytical expressions for ⟨p 2 (t)⟩ and ⟨x 2 (t)⟩.Performing the Laplace transformation of Caputo derivative for n In order to calculate the solution of p(t) and x(t) of the heavy particle for a super-diffusive process, i.e., when 1 < α ≤ 2 and 1 < β ≤ 2, one can take the Laplace transform Eqs. ( 2), (3) and using Eq. ( 9), it can be obtained easily as Now, substituting Eq. ( 11) into Eq.( 10), for simplicity, we take M = 1 [94], we obtain Taking the inverse Laplace transform of Eqs. ( 11) and ( 12), we obtain and where is the two-parameter Mittag-Leffler function for α > 0, β > 0, and z is a complex number [95] and the Laplace transform of the two-parameter Mittag-Leffler function is [96], Next we calculate ⟨p 2 (t)⟩, and ⟨x 2 (t)⟩ analytically for the Brownian particle.

Purely diffusive motion
We are calculating the ⟨p 2 (t)⟩, and ⟨x 2 (t)⟩ for a purely diffusive (when γ = 0 in the FLE), 1-D motion.Then Eq. (3) simplified to with initial momentum, p(t 0 ) = 0. Then to calculate ⟨p 2 (t)⟩ we follows Eq. ( 13) and Eq. ( 16), then we get Solving further Eq.( 17), using the correlation properties defined in Eq. ( 4), we have, Similarly, one can calculate analytical solutions of ⟨x 2 (t)⟩ for γ = 0, and with initial conditions, x(t 0 ) = 0, using Eqs.( 14), ( 16) as We have plotted Fig. 1 for the case of purely diffusive motion, we depict the variation of ⟨p 2 (t)⟩ over time (left panel) and ⟨x 2 (t)⟩ (right panel), each for different values of α and β, while keeping D constant at 0.1 GeV 2 /fm.For simplicity, a non-relativistic case has been considered with M = 1 GeV.In Fig. 1 (left panel), results obtained from Eq. (19) show that ⟨p 2 (t)⟩ exhibits nearly linear evolution with time when α → 1 (as shown by black solid lines), which is giving ⟨p 2 (t)⟩ = 2Dt.In this scenario, the Caputo fractional derivative reduces to an ordinary first-order derivative; in this limit of α, the FLE is reduced back to the LE.On the other hand, when the value of α > 1, the linear increase in ⟨p 2 (t)⟩ with time converts into nonlinear growth over time.It becomes evident that the diffusion process evolves more gradually during the initial time, particularly up to 1.5 fm/c.In contrast, around 2 fm/c, the pattern is reversed, now a higher value of α, giving a faster diffusion.Similarly, with the same input parameter, Fig. 1 (right panel) results obtained from Eq. ( 20) depicts the mean square displacement, ⟨x 2 (t)⟩ varies with time.In the limit, α, β → 1(as shown by black solid lines), ⟨x 2 (t)⟩ varies with t 3 (as shown by black solid lines), for the normal diffusion.In this case, we considered the scenario of purely diffusive motion (γ = 0 in the FLE), where the diffusion of the Brownian particle predominates due to this approximation.Consequently, ⟨x(t) 2 ⟩ of the non-relativistic Brownian particle varies with t 3 for pure normal diffusion.It is crucial to emphasize that because of γ= 0, the variation of ⟨x(t) 2 ⟩ with t 3 is unrelated to Eq (1).Given the limiting nature of this scenario, the calculated results seem to contradict the expected behaviour of normal diffusion, as expected in Eq. ( 1).Once our analytical and numerical calculations align, we will incorporate the drag term into the relativistic Langevin equation for the charm quark in the subsequent subsection.The reason for choosing this limiting case will be discussed in the following subsection.With the higher values of α and β, ⟨x 2 (t)⟩ varies with the power of time greater than cubic power.

B. A numerical scheme for the fractional Langevin equation
In classical stochastic differential equations driven by Brownian motion, the Itô or Stratonovich stochastic calculus is typically employed for solutions.However, these methods are not applicable to the FLE driven by fractional Brownian motion, as it is not a semimartingale (see Ref. [97] for detailed proof).Although the Monte Carlo method is a reliable solution for stochastic differential equations, it is not appropriate in this case because it relies on independent sequences.Still, the sequences in fractional noise are dependent.Therefore, we have employed L2 numerical schemes, as detailed in [98].These schemes are among the most effective for discretizing the Caputo fractional derivative and provide a key contribution to this article.We numerically solve the FLE in Eqs.(2), (3).To validate our numerical computations, we compared the numerical results with the analytical ones for the 1-D FLE, considering γ = 0.This step is crucial to ensure the accuracy of our numerical method, which will be employed in solving the relativistic Langevin equation that will be discussed in the following subsection.To solve the FLE, numerical methods are commonly catego-rized as indirect or direct.Since time-fractional differential equations can generally be reformulated into integrodifferential equations, solving such equations corresponds to indirect methods.On the other hand, direct methods focus on approximating the time-fractional derivative itself.This aspect constitutes the core of our work, where our goal is to discretize the fractional derivative directly without transforming the associated differential equation into its integral form.Our numerical algorithm for solving the FLE is based on the three-step scheme using the central difference method shown in Appendix VI.

The fractional Langevin equation of relativistic
Brownian particles: QGP We extended the solution of FLE for the nonrelativistic Brownian pariticle as defined in Eqs. ( 2),(3), to describe the HQs dynamics in the QGP medium in the relativistic limit as follows, where, E = p 2 + M 2 represents the energy, and p denotes the momentum of the HQs.The γ can be related to the diffusion coefficient through an FDT as [99][100][101], Since both Eq. ( 21) and Eq. ( 22) are non-linear differential integral equations, preventing analytical solutions through methods such as Laplace transformations, as done for the non-relativistic case for Eqs. ( 2), (3).The solution of relativistic FLE defined in Eq. ( 22) and Eq. ( 21) is possible solely through numerical methods.We utilize the numerical approach to compute ⟨x 2 (t)⟩, ⟨p 2 (t)⟩, the momentum distribution, and R AA of the HQs.For superdiffusion, the relativistic FLE can be written in the discrete form using Eq. ( 30) ( given in appendix VI) as follows, and where Here t 0 is the initial time, and ∆t is the time step.

A. A check for non-relativistic
To analyze the numerical scheme of FLE, we compute ⟨x 2 (t)⟩ from Eq. ( 24) and ⟨p 2 (t)⟩ from Eq. ( 25), subsequently comparing the results with the analytical solutions derived from Eq. ( 20) and Eq. ( 19), respectively.The FLE has been solved for the 1-D motion of the heavy particle within the nonrelativistic limit, as explained in Sec.II.For illustrative purposes, we have considered constant diffusion coefficient, D = 0.1 GeV 2 /fm and mass of the heavy particle, M = 1 GeV and γ = 0 (i.e., pure diffusion case) in FLE.In the left panel of Fig. 1, the numerically computed ⟨p 2 (t)⟩ for α = 1.001 (brown dashed line), α = 1.2 (grey dashed line), α = 1.4 (red dashed line), and α = 1.6 (orange dashed line) is presented.Notably, these numerical results align with the analytic solutions from Eq. ( 19) and match the other values of the α.For an additional test of our numerical approach, we have calculated ⟨x 2 (t)⟩ of the heavy particle as shown in Fig. 1 (right panel).In the same figure, the numerically calculated ⟨x 2 (t)⟩ for α = β = 1.001 (brown dashed line) and α = β = 1.2 (grey dashed line), α = β = 1.4 (red dashed line), and α = β = 1.6 (orange dashed line) is displayed.Again, the numerical results match the analytic results from Eq. (20).One can notice that the numerical simulation agrees with the analytical result, indicating that our numerical scheme works correctly.
B. The evolution of ⟨p 2 (t)⟩ and ⟨x 2 (t)⟩ of the HQs.
For our analysis, the definition we used for ⟨p 2 (t)⟩ is given by We have computed ⟨p 2 (t)⟩ over time for four different values of α, namely α = 1.001 (black line), α = 1.2 (red line), α = 1.4 (green line), and α = 1.6 (blue line) at T = 250 MeV, M = 1.3 GeV and D = 0.1 GeV 2 /fm correspond to the charm quark.The corresponding results are shown in Fig. 2. In this scenario, the initial momentum is set to be p x (t 0 ) = p y (t 0 ) = 0.
From the left panel of Fig. 2, it is evident that as the magnitude of α increases, the behaviour of the process indicates superdiffusion for the charm quark within the hot QCD medium.The impact of superdiffusion becomes more noticeable at higher values of α.It is noticeable that as α → 1, the superdiffusion process converges back to normal diffusion.Additionally, in the later stages, the mean squared momentum ⟨p 2 (t)⟩ tends to approach 3M T [100], and also the FLE defined in Eq. ( 25) simplifies to the standard LE described in Ref. [39,40,43,45,100].This transition emphasizes the connection between superdiffusive behaviour governed by the FLE and the conventional diffusion process characterized by the standard LE.In the right panel of Fig. 2, we present a subset of the results depicted in the left panel, focusing on the earlytime evolution of ⟨p 2 (t)⟩ for four different values of α.It is worth noting that in the initial time, the diffusion process gradually evolves for higher values of α.However, for times beyond 3 fm/c, the trend reverses, and larger α corresponds to a faster diffusion, as anticipated in Fig. 2 (left panel).In Fig. 3, we have calculated the evolution of ⟨x 2 (t)⟩ of the relativistic charm quark over time for three values of α = β = 1.001, 1.2, 1.4, maintaining other parameters consistent with those illustrated in Fig. 2. The initial conditions are set to be x(t 0 ) = y(t 0 ) = 0.It can be noticed that both α = β = 1.2, 1.4, a distinct shift towards superdiffusion is observed (shown in Fig. 3) as we discussed initially in Eq. (1).The larger values of α and β contribute to this notable change in behaviour, underlining the complex dynamics of the system.As α → 1 and β → 1, the behaviour reflects normal diffusion.In this limit, at a later time, ⟨x 2 (t)⟩ exhibits a proportional relationship with t (same as in Eq. ( 1)), as described in [100,102].This puzzling behaviour of ⟨p 2 (t)⟩ and ⟨x 2 (t)⟩ due to the anomalous diffusion, which was not explained before in the context of the HQ dynamics in a hot QCD medium.
In the following sections, we have calculated R AA and the charm quark momentum distribution, dN/dp T , to see the effect of superdiffusion processes.To determine the interaction of the HQ with the thermalized bath consisting of massless quarks and gluons, we employ perturbative Quantum Chromodynamics (pQCD) transport coefficients for elastic processes with the wellestablished diffusion coefficients [102].Here, the Debye mass, m D = g T T screens the infrared divergence associated with the t-channel diagrams.

C. Nuclear modification factor
To analyze the impact of the superdiffusion on the experimental observable, we have calculated the nuclear suppression factor, R AA (p T ) for the HQs, which is defined as follows [100], The momentum spectrum, f τ f (p), of charm quarks calculated for the time evolution, τ f = 6 fm/c in our computational results and f τi (p) is initial momentum distribution of the charm quark.The initial momentum spectra, f τi (p), is taken according to the fixed order + next-toleading log (FONLL) calculations, which has been shown to be capable of reproducing the spectra of D-mesons produced in proton-proton collisions through fragmentation [103,104], the initial momentum spectrum is written as, where the parameters are estimated as follow; x 0 = 6.365480 × 10 8 , x 1 = 9.0 and x 2 = 10.27890.R AA (p T ) ̸ = 1 implies that charm quarks undergo interactions with the medium.These interactions lead to modifications in the spectrum of charm quarks.In Fig. 4, the behaviour of R AA is depicted as a function of p T , which are calculated using the FLE as defined in Eq. ( 25) for different values of α.The calculations are performed at two distinct temperatures, T = 250 MeV (left panel) and T = 350 MeV (right panel).At T = 250 MeV, the dominating influence is the drag force across the entire range of p T .As p T increases, the significance of energy loss becomes more pronounced, which can be noticed in Fig. 4 (left panel).When α > 1, the normal diffusion converts into superdiffusion, it is evident for α = 1.2 (red line), α = 1.4 (blue line), and α = 1.6 (green line).With an increasing value of α, there is a notable decrease in magnitude R AA (more suppression) at high p T .For α → 1 (depicted by the black line), the behaviour of the R AA aligns with normal diffusion, consistent with findings available in the literature for the same input parameters [40,43,46].Fig. 4 (right panel), corresponds to T = 350 MeV, the observed behaviour is attributed to the diffusion-dominated propagation of the HQs within the hot QCD medium.This dominance of diffusion mechanisms effectively leads to the diffusion of low-momentum charm quarks to higher momentum states.The high temperature, T = 350 MeV, enhances the significance of diffusion processes, resulting in a distinct pattern in the R AA as compared to R AA at T = 250 MeV.For α > 1, a significant reduction in the magnitude of R AA is observed, indicating more suppression at high p T .Notably, for the highest considered value of α = 1.6, a larger proportion of particles tends to remain at low p T , a consequence of the superdiffusion process.This behaviour underscores the complex dynamics associated with superdiffusion and its impact on the R AA of the HQs in the QGP medium.

D. Momentum spread of HQs
We show the evolution of charm momentum distribution, dN/dp T , at static temperatures, T = 250 MeV (top panel) and T = 350 MeV (bottom panel) in Fig. 5.The dN/dp T evolution is performed for various values of α at a final evolution time of τ f = 6 fm/c.To understand the impact of superdiffusion on charm quarks within the QGP, we take initial conditions where all charm quarks are concentrated within an extremely narrow p T bin, creating a delta-like distribution at p(t 0 ) = 10 GeV (magenta line).It is observed that the interaction of the HQs with the QGP medium results in the spreading of dN/dp T .Subsequently, the evolution of this distribution is analyzed using the FLE as defined in Eq. (25).We have observed the evolution of dN/dp T at higher values of α for the case of superdiffusion, where, α = 1.2 (red line), α = 1.4 (blue line), and α = 1.6 (green line).As depicted in Fig. 5, when α > 1, the distribution dN/dp T undergoes a notable spread and average momentum shifts towards lower values of p T under the influence of diffusion and drag coefficients, respectively.Specifically, for α = 1.6 (represented by the blue line), the extent of spreading is more pronounced, and the average momentum shifts towards the lower p T compared to other α values such as 1.001, 1.2, and 1.4.However, for α = 1.001 (black line), corresponding to the normal diffusion coefficient as explained in Ref. [41,43].At the same time, the total area under the curve remains constant for all values of α

IV. CONCLUSION AND OUTLOOKS
In this paper, we have discussed anomalous diffusion through the FLE with the Caputo fractional derivative, specifically focusing on superdiffusion in the context of the HQs in the QGP medium.Notably, the mean squared displacement of the particle exhibits a power-law dependence on time (as shown in Fig. 3).Our analysis discussed the scenarios where the values of α and β were taken as 1.001, 1.2, 1.4, and 1.6, showcasing the effects of superdiffusion.We have demonstrated that as α approaches 1, the superdiffusion reverts to normal diffusion, verified by numerical and analytical calculations for ⟨x 2 (t)⟩ and ⟨p 2 (t)⟩ in the non-relativistic 1-D case (see Fig. 1).Initially, in Section II.A.1, we focused on purely diffusive motion with γ = 0 in the FLE.Because of the dominance of the diffusion term, ⟨x(t) 2 ⟩ varies with t 3 , contrary to normal diffusion expectations.Upon aligning analytical and numerical results, then we later incorporated γ into relativistic FLE, confirming normal diffu-sion for relativistic charm quarks (as discussed in Eq. ( 1)).Several key quantities characterizing the dynamics of the HQs under superdiffusion have been computed, including the ⟨p 2 (t)⟩, ⟨x 2 (t)⟩.Extending our analysis, we incorporated physical observables R AA of the HQs in the QGP medium.We then shifted our focus to the momentum spread, dN/dp T , utilizing an initial momentum distribution at p(t 0 ) = 10 GeV.The FLE, with the HQs moving under dissipative and random forces, was solved with transport coefficients serving as input parameters.Our findings indicate that superdiffusion results in more suppression in R AA .To select specific values for α and β within the ranges of 1 < α ≤ 2 and 1 < β ≤ 2, it is essential to be able to simultaneously describe the experimental observables, such as R AA and v 2 , for the entire measured range of p T .Initially, the formation of a small R AA (signifying strong suppression) can occur rapidly at the beginning of the QGP, while the development of significant v 2 is more sensitive to later stages of evolution.Consequently, substantial interactions may not coincide with a significant build-up of v 2 since the bulk medium has yet to establish significant elliptic flow.Exploring anomalous diffusion may offer the potential to generate notable v 2 magnitude, possibly leading to intensified interactions between the HQs and the evolving bulk at later stages.However, a thorough analysis of v 2 in the presence of anomalous diffusion requires a refined study incorporating realistic initial conditions, including the initial geometry and expansion of the fireball.Subsequent comparison with experimental data will be crucial for validating and refining these findings, representing an essential aspect of future investigations in this study.
In the future, we plan to explore the impact of superdiffusion on the HQ dynamics in the QGP medium, especially considering the time correlation of thermal noise.This study devotes the groundwork for more realistic conditions, which should incorporate an exact initial geometry and an expanding medium in the near future.Superdiffusion might impact various observables, such as HQs directed flow, particle correlations, etc.Given the simplifying assumptions made in our current study, it is challenging to anticipate the specific modifications.A more comprehensive and quantitative analysis will be performed in near future investigations to understand these phenomena further.

V. ACKNOWLEDGMENTS
I gratefully acknowledge Dr. Santosh Kumar Das for his invaluable advice, Aditi Tomar for insightful discussions and inspirational encouragement, and Mohammad Yousuf Jamal for his numerous informative contributions, collectively enriching the content of this paper.Additionally, I acknowledge the partial support from DAE-BRNS, India, Grant No. 57/14/02/ 2021-BRNS.I also acknowledge partial support from the SERB Fellowship Project Code No.RD/0122-SERBF30-001.

VI. APPENDIX
Consider a partition {t n = nT N , ; 0 ≤ n ≤ N } of time interval [0, T ].The three-step numerical scheme called L2 approximation is used for the case of superdiffusion.The second derivative u (2) is approximated using a central difference formula, and the resulting numerical scheme involves the values of u at the previous three-time points u n−1 , u n−2 , and u n−3 .29) where the coefficients b j = (j+1) 2−ν −j 2−ν Γ(3−ν)∆t ν in the scheme is determined by the difference formula for the second derivative and are used to account for the fractional order.
FIG. 1: ⟨p 2 (t)⟩ versus time (left panel) and ⟨x 2 (t)⟩ versus time (right panel) for a 1-D, purely diffusive motion for different values of α and β.Analytic solutions are represented by solid lines, and numerical (Numer.)solutions are represented by dashed lines.

6 FIG. 4 :
FIG. 4: The R AA as a function of p T , with t = 6 fm/c, and at two different temperatures (T = 250 MeV (left panel) and T = 350 MeV (right panel) with four different values of α.

FIG. 5 :
FIG. 5: The evolution of charm quark momentum distribution at τ f = 6 fm/c as a function of p T , for various α and at two different temperatures (T = 250 MeV (above panel) and T = 350 MeV (below panel).Assuming an initial delta distribution centered at p(t 0 ) = 10 GeV.