Note that in Eq. (2.15) the memory integral diverges for H smaller than 1 / 2, such that the Hurst exponent in the FLE is restricted to the range 1 / 2 < H < 1. It can be shown that the relaxation dynamics governed by the FLE (2.15) follows the form 4 generated by the Hosking method which is known to be an exact but time-consuming algorithm [56]. We checked that in the one-dimensional case the generated FBM in free space successfully reproduces the theoretically ex- pected behavior, the mean squared displacement (1.1), the fractal dimension d$_{f}$ = 2 − α/ 2 of the resulting trajec- tory, and the first passage time distribution. To simulate the confined motion, reflecting walls were considered at locations ± L for each coordinate. For instance in the 1D case, if | x $^{H}$( t ) | > L at some time t , the particle bounces back to the position x $^{H}$( t ) − 2 | x $^{H}$( t ) − sign( x $^{H}$) L | . Sim- ilar reflecting conditions were taken into account in the multi-dimensional case. for the first moment, where v$_{0}$ is the initial particle veloc- ity. The rescaled friction coefficient is γ = γ Γ(2 H − 1) /m . The coordinate variance behaves as where 〈 v 2 $_{0}$〉 = k$_{B}$T/m is assumed, and we employed the generalized Mittag-Leffler function [53] whose asymptotic behavior for large z is Thus the mean squared displacement shows a turnover from short time ballistic motion to long time anomalous diffusion of the form [54, 55] Therefore, for persistent noise with 1 / 2 < H < 1 the re- sulting motion is in fact subdiffusive, i.e., the persistence of the noise has the opposite effect than in FBM. In analogy to our discussion of FBM in a d -dimensional embedding the FLE is generalized to where y ( t ) = ∑ $_{i}$y$_{i}$ ( t )ˆ x$_{i}$ , ξ ( t ) = ∑ $_{i}$ξ H i ( t )ˆ x$_{i}$ , and ⃗ K is the memory tensor which is in diagonal form (i.e., K ij = K ( t − t $^{′}$) δ$_{ij}$ ) in the absence of motional coupling between different coordinates. We here briefly review the simulations scheme used to produce time series for FBM and FLE motion. d -dimensional FBM is simulated via Eq. (2.7) by nu- merical integration of ξ H i ( t ). The underlying FGN was where dh is the time increment. When evaluating Eq. (3.2), a reflecting boundary condition was considered in the sense that y$_{n}$ → y$_{n}$ − 2 | y$_{n}$ − sign( y$_{n}$ ) L | if | y$_{n}$ | > L . To show the reliability of our simulation, we compare the simulation result with the well-known solution for free space motion. In Fig. 1, we simulate the subdiffu- sion case with the parameters values, H = 5 / 8, m = 1, v$_{0}$ = 1, y$_{0}$ = 0, γ = 10, k$_{B}$T = 1, and dh = 0 . 01. From 200 simulated trajectories, we obtain the ensemble av- eraged and time averaged mean squared displacements, 〈 y $^{2}$( t ) 〉 and δ $^{2}$(Δ , T ), and compare them with the exact solution Eq. (2.19). Note that 〈 y $^{2}$( t ) 〉 should be identical to δ $^{2}$(Δ , T ), Eq. (5.2), with t regarded as the lag time Δ due to the ergodicity of the FLE motion in free space [25]. The deviation from the exact solution is markedly re- duced with decreasing time increment dh . With our cho- sen value dh = 0 . 01 the mean squared displacements ob- tained from simulation appear to be in good agreement with the theory. where v ( t = 0) = v$_{0}$ . This stochastic integral equa- tion can be evaluated by the predictor-corrector algo- rithm presented in Ref. [57] with the FBM x $^{H}$( t ) inde- pendently obtained by the Hosking method. We calcu- lated y ( t ) = y$_{0}$ + ∫ t 0 v ( t $^{′}$) dt ′ by the trapezoidal algorithm. For discrete time steps, the equation of motion is given by In simulating FLE motion, we follow the numerical method presented by Deng and Barkai [25]. First, inte- grating Eq. (2.15) from 0 to t , we obtain the Volterra integral equation for velocity field v ( t ) = dy ( t ) /dt B. Fractional Langevin equation motion v ( t ) = − γ (2 H − 1) m ∫ t 0 ( t − t $^{′}$) 2 H − $^{1}$v ( t $^{′}$) dt ′ + v$_{0}$ + η m x $^{H}$( t ) , (3.1) y$_{n}$$_{+1}$ = y$_{0}$ + dh 2 ( v$_{0}$ + v$_{n}$$_{+1}$ ) + dh n ∑ i =1 v$_{i}$, = dh 2 ( v$_{n}$ + v$_{n}$$_{+1}$ ) + y$_{n}$, (3.2) m d 2 y ( t ) dt 2 = − γ ∫ t 0 ⃗ K · d y dt $_{′}$dt ′ + ξ ( t ) , (2.23) 〈 y $^{2}$( t ) 〉 ∼ { t $^{2}$, t → 0 t 2 − 2 $^{H}$, t → ∞ . (2.22) E$_{α,β}$ ( z ) = − ∞ ∑ n =1 z − n Γ( β − αn ) . (2.21) E$_{α,β}$ ( z ) = ∞ ∑ n =0 z n Γ( αn + β ) (2.20) 〈 y $^{2}$( t ) 〉 = 2 k$_{B}$T m t $^{2}$E$_{2}$$_{H,}$$_{3}$ ( − γt 2 H ) (2.19) 〈 y ( t ) 〉 = v$_{0}$tE$_{2}$$_{H,}$$_{2}$ ( − γt 2 H ) (2.18) III. SIMULATIONS SCHEME A. Fractional Brownian motion