Subdiffusion in an array of solid obstacles

More than a decade ago, Goychuk reported on a universal behavior of subdiffusive motion (as described by the generalized Langevin equation) in a one-dimensional bounded periodic potential (Goychuk 2009 Phys. Rev. E 80 046125) where the numerical findings show that the long-time behavior of the mean squared displacement is not influenced by the potential, so that the behavior in the potential, under homogenization, is the same as in its absence. This property may break down if the potential is unbounded. In the present work, we report on the results of simulations of subdiffusion in a two-dimensional (2D) periodic array of solid obstacles (i.e. in an unbounded potential) with different packing fractions. It is revealed that the universal subdiffusive behavior at long times is not influenced by the presence of solid scatterers, whose presence influences the behavior at intermediate times only. This result is discussed as having possible relations to the emerging problem of interpretation of results on trajectories of tracers spreading in the brain’s extracellular space.

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The motion of microscopic tracers in crowded liquids (artificial media or interiors of biological cells) and complex-structured soft matter (e.g.hydrogels) often shows anomalous diffusion [1,2].Such behavior can originate from different physical mechanisms related both to the structure of the inhomogeneous medium in which the tracer's walk occurs, and to the properties of a carrier fluid filling this structure [3].Among this variety, we will focus in this work on a particular example of viscoelastic medium that additionally contains obstacles.Recently, this case gained a special interest due to the modern development of experimental tools based on single-particle tracking, which allows obtaining detailed data for simultaneous quantification of transport and biomechanical properties of living matter [4,5].In many cases, generalized Langevin equations deliver an instrument for the adequate description of the corresponding phenomena.
It should be pointed out that the majority of the corresponding works deals with the exploration of the cellular interior.At the same time, there exists another realm of biophysical phenomena, which is significantly less explored: the spread of high-molecular-weight markers and nanoparticles in the brain's parenchyma.There is an emerging understanding that the corresponding processes can include not only normal Fickian (Brownian) diffusion (sometimes with a position-dependent diffusion coefficient) but also Brownian-yet-not-Gaussian diffusion, and subdiffusion, see e.g.[6] for a review.
The subdiffusive spread is especially typical for the case of nanoparticles [7,8] which have relatively large sizes comparable with gaps forming the extracellular space (ECS), and includes effects of trapping, temporal binding, etc leading to a hindered motion.It is worth noting, however, that the brain's ECS is not just a system of interconnected disordered voids filled with the liquid with properties similar to the cerebrospinal fluid.The parenchyma contains also a fine macromolecular network composed of fibrous proteins and polysaccharides [9][10][11].Such a composition of the carrier medium is well-known as affecting the spread of particles [12] and leads to viscoelastic properties similar to those of polymer solutions [13].The interaction of walking nanoparticles with such medium can effectively activate viscoelastic retarding, e.g.via a local swelling, which is crucial for understanding the controlled drug delivery into the brain's parenchyma [14][15][16].
It should be pointed out that when considering the spread of macromolecular or nanoparticle tracers in the brain parenchyma, one should take into account an interplay between microscopic and macroscopic scales [17][18][19].The conventionally used markers cannot penetrate the walls of surrounding cells, i.e. one needs to address the microscopic problem of a random walk in a system of reflecting obstacles within the borders of a complex ECS.At the same time, the widely used biophysical methods aimed at the characterization of the ECS state often operate at the macroscopic level, especially in vivo, i.e. they deal with either the spatiotemporal evolution of the markers' cloud or with the effective diffusion coefficient determined on the base of the former.This is known as a problem of the interplay between the porosity (packing fraction of hard obstacles) and tortuosity (average length of the paths between the source and the registration site).Both lead to lowering the effective diffusion coefficient compared to the one in the absence of obstacles.The corresponding dependence is not trivial even for the case of the normal Brownian motion [20,21]; in some experiments, however, subdiffusion is observed [14,22,23], but this situation is still underexplored despite the emergence of new results of single-particle tracking in the brain's ECS [24].However, such studies are still scarce.
Another recently emerging field of study related to the interplay of the brain's anatomic features and models of fractional random walks in a medium with obstacles, unrelated to diffusive spread, is the analogy between the shape of serotogenic axons and the random walker's trajectories.To simulate the shape and local density of axons, the authors of [25,26] modeled them by trajectories of the fractional Brownian motion (fBm) reflected by impenetrable curved boundaries and small, not very densely packed, scatterers.It should be pointed out that such a discussion poses a specific problem of adequate simulation of reflections at the boundaries [27].The approaches when one simply rejects steps leading to the impenetrable region or introduces artificial, not a perfectly elastic walls are suitable only for the case of large empty reservoirs or sparsely located obstacles.The opposite case of densely packed obstacles is more complicated and, to the best of our knowledge, unexplored up to date.The problem of implementation of hard boundaries persists also for the generalized Langevin equation approaches which may show a memory structure different from those of fBm.
Motivated by a class of problems discussed above, in the present work we consider a subdiffusive motion of particles as described by a generalized Langevin equation with a power-law memory kernel in an array of impenetrable solid obstacles.The aim is to identify the dependence of the effective subdiffusion's parameters which can be accessed from both macro-and microscopic measurements (like the exponent of the subdiffusion and the subdiffusion coefficient) on the packing fraction of the obstacles.Note that such a problem implies the consideration of an at least a 2D case, in which the walker has the possibility to path between the obstancles and to escape to infinity.In one dimension, impenetrable walls will always lead to confinement of the walker [28,29], which corresponds to a very different situation.Our discussion here concentrates on the simplest situation corresponding to homogenization of subdiffusion in a periodic array of obstacles.Although our approach can be also used for other arrangements, the periodic array provides a convenient picture for both numerical simulations and for the discussion of the interplay between the cases of unbounded and bounded potentials.Thus, it is the first step on the way to understanding more complex, biologically relevant, situations, showing some important effects in the simplest setup and in their absolute purity.Our main interest is devoted to what happens at relatively high packing fractions which are however far enough from the cases when narrow escape properties start to play a role.The implementation of the generalized Langevin's simulations in the case of subdiffusion in presence of multiple closely located impenetrable walls is a numerical task which is usually computationally complicated, but which can be effectively solved by a combination of Markovian embedding as proposed by Goychuk, and modeling of impenetrable obstacles with 'a background field' trick which allows for easy implementation of specular reflection from curved surfaces.The main result of our simulations is that the presence of the obstacles does not change the long-time properties of the subdiffusion.In this sense our 2D situation is similar to the one for one-dimensional motion in a bounded periodic potential as considered in [30].

The generalized Langevin equation in equilibrium bath
The generalized Langevin equation (GLE) for the motion of a particle in the presence of an equilibrium bath reads: with M being the mass of the particle, K(t) the memory kernel, f(x) the external force assumed conservative with potential U(x), f(x) = −∇U(x), and ζ(t) being the noise term.The noise is assumed to be Gaussian, and to have zero mean.The coordinate x of a particle may be considered as a scalar in a one-dimensional problem or as a vector in higher dimensions (the same is true for the force f and the noise (random force) ζ).In general, in more than one dimension the memory kernel K(τ ) is a tensor-valued function of time; in what follows we will consider only the isotropic situation in which this is represented by a diagonal matrix with equal diagonal elements K(τ ).
The equation represents quite a general form of interaction of a tagged particle with the bath, can be derived from the underlying Hamiltonian dynamics with the help of projection operator formalism (see e.g.[31]), and can also be derived explicitly for specific bath models (e.g. for Kac-Zwanzig bath, see [32]).The property of the bath to be at equilibrium at temperature T leads to the fluctuation-dissipation relation between the memory kernel and the noise correlation function, (here in a one-dimensional notation) for τ ⩾ 0, where the average can be considered both as a time and as an ensemble average due to the assumed ergodicity of the noise.In the isotropic 2D situation considered in section 4 ff this relation holds for each Cartesian component separately.Different Cartesian components of the noise are uncorrelated.For τ < 0 the memory kernel K(τ ) vanishes identically, while is an even function of τ .The kernel K(τ ) will be assumed to be a continuous function of τ for all cases except for the diffusive one (with the white noise), when it tends to a δ-function.
Depending on the exact form of the kernel, the equation can describe a variety of diffusive motions, but also sub-and superdiffusion, and is of wide use, see e.g.[28] and references therein.In the present work we concentrate on the subdiffusive case.

The force-free case
Before passing to the general situation, let us discuss the motion of the particle in one dimension in the absence of the external potential.Our notation and discussion here follow to a large extent the ones in [32].The detailed presentation of known results will be of use when discussing the role of external force, and the parameters used in numerical simulations of the 2D situation which will be the main topic of our work.The generalized Langevin equation in one dimension reads To analyze the diffusion behavior we will consider this equation with initial conditions x(0) = x 0 , ẋ(0) = v 0 , where v 0 will later be assumed to be a random variable following the Maxwell-Boltzmann distribution at temperature T.
Transforming equation (3) to the Laplace domain with and readily obtains the solution Defining the Green's function of the GLE H(t) with the Laplace representation one may write the solution in the time domain as We moreover note that according to the initial value theorem for all cases when K(s) does not tend to −Ms for s → ∞.Since for t → 0 K(t) → ⟨ζ 2 (t)⟩ > 0 and K(t) in our model is assumed continuous (except for the limit of the white noise), its Laplace transform for large s is non-negative, and the situation cannot be realized.Now we start with x 0 = 0 and consider the mean squared displacement (MSD) of the particle from its initial position where the function H(t) is non-random, and the averages are taken over the realizations of the noise and over the initial conditions for v 0 .If we start from the thermalized conditions ⟨Mv 2 0 ⟩ = k B T (here in one dimension) then the first term equals to Mk B TH 2 (t).If we start with zero velocity, it vanishes.The second term vanishes identically when we start at zero velocity, and vanishes on the average in the thermalized case if v 0 and ζ(t) are independent (i.e. when the particle and the bath are independently set to equilibrium at t = 0, as we will do in our simulations), and ⟨v 0 ⟩ = 0.In all other cases it gives a decaying, subleading contribution.Rewriting the third term as and using equation ( 2) and the symmetry of the correlation function of the noise we further transform the expression: The internal integral G(ξ) = ´ξ 0 K(ξ − η)H(η)dη is a convolution, and can be calculated in the Laplace domain: In time domain this means where the dot denotes the time derivative.Now the whole integral reads Combining this with the first term in equation ( 5) we see that In the multidimensional case with independent noises in different directions this expression holds for each direction separately.Now let us assume that R 2 (t) ∝ t α , namely that for t sufficiently long.The case 0 < α < 1 corresponds to subdiffusion, the case α = 1 to normal diffusion, and the case 1 < α ⩽ 2 to superdiffusion (as we proceed to show, the GLE with equilibrium bath does not describe generic superballistic regimes).Using equation ( 8) we now obtain that for long times and therefore for s sufficiently small.Assuming the behavior of H(s) as given by equation ( 4), we see that and therefore K(s) represents the leading term for all 0 < α < 2. For 0 < α < 1 one has K(s) ∼ s α−1 , and the integral K(s = 0) = ´∞ 0 K(t)dt diverges.For α = 1 the value of K(s) tends to a constant, and the corresponding integral converges to a finite value.
In the case of subdiffusion, the immediate application of the Tauberian theorem gives us the behavior of K(t) in the time domain: for long times one has For α = 1 the kernel can be approximated by a δ-function: In the case of superdiffusion 1 < α < 2 one has K(s) ∼ s α−1 → 0, and the integral K(s = 0) = ´∞ 0 K(t)dt vanishes.In this case one has to consider the integral I(t) = ´t 0 K(t ′ )dt ′ with I(∞) = 0. Since K(0) > 0, this may only be the case if the function K(t) changes its sign at least once.This can be also seen directly from equation (11), which, as one can show, is still applicable for this case, and gives negative values of K(t) for this case for t long.The fact that a particular integral has to vanish identically in the long time limit makes long-time subdiffusion fragile with respect to external perturbations under which it easily changes to a long-time diffusive behavior.
In the case of superballistic behavior, α > 2, equation (10) gives In what follows we concentrate on the case of subdisffusion 0 < α < 1 and note that the motion with the subdiffusion coefficient D * corresponds to the following relation between D * , defining the long-time MSD, equation ( 9) and K 0 being the prefactor of the power-law in

(Non)-influence of the external potential
In [30], Goychuk made a remarkable finding that the long-time subdiffusive behavior of a particle in a bounded periodic potential is not influenced by the strength of this potential.The potential strength only governs the crossover time to the long-time subdiffusive asymptotics, being the same for all potential strengths (and the same as without the external potential, with R 2 (t) ≃ D * t α ).In other words, this means that the a subdiffusion in a periodic potential is, under homogenization, exactly the same as if no potential was acting on the trancer particle.On the other hand, in a random Gaussian potential [33] with differently parameterized spatial correlation functions (e.g.decaying exponentially or according to a power-law as representative for structures of some polymer and biophysical disordered media, see [34,35] and references therein) the presence of such a potential may lead to the renormalization of the subdiffusion coefficient D * as well as to an extremely rich intermediate-time behavior.There is no regular theory describing the effects, but a posteriory one can figure out, when the deviations from the free behavior may appear.One of them corresponds to unbounded potentials.The idea behind our qualitative discussion (lacking mathematical rigor) is connected with fast local equilibration of the PDF in a bounded potential, so that the PDF of the particles' positions at time t in such a subdiffusion is well-described by the Gaussian modulated by the Boltzmann factor: with N being the correction to the normalization factor corresponding to a Gaussian envelope, and D * is an effective subdiffusion coefficient, which (in principle) may differ from D * .The idea behind this representation is the same as for non-confining potentials in the Markov case, see [36].Moreover, one has to bear in mind the fact the effective exponent α also may in principle differ from the one in the free motion, but here we are considering the situations (diffusion and subdiffusion) when it does not.The corresponding PDF for a sinusoidal force f(x) = −f 0 sin 2π x L with period L, i.e. in a cosinusoidal potential as obtained in simulations, is shown in figure 1.These simulations present a numerical evidence for equation ( 13), but equation ( 13) can also be made plausible by considering relaxation timescales like it is done in [30].The simulations show that the long-time form, equation (13), is achieved independently from the initial position of the tracer in the first cell and its initial velocity, i.e. that the relaxation of the particle within the cell is essentially a fast process compared to the spread.One also can note that the MSD at very long times given by the distribution equation ( 13) will tend to ⟨(x − x 0 ) 2 ⟩ = D * t α defining the envelope's width.This can be seen by splitting the integration domain into intervals corresponding to the periods of potential, applying the mean value theorem for integrals to each of them and carefully considering the ensuing sum over the points sampling a Gaussian, one per interval in the case when the envelope function within a single interval hardly changes.
Note that in the course of the time, for D * t α ≫ L, the approximation, equation ( 13), would give us the stationary distribution of f(t) = f[x(t)] which will simply correspond to its distribution over the single period of the potential in thermodynamic equilibrium, so that the force f (x) essentially reaches equilibrium in spite of the fact that x does not have an equilibrium distribution.This, and several other properties of the situation are explicitly discussed in the appendix for the special case of a cosinusoidal potential, equation (14), corresponding to the sinusoidal force f(x) = −f 0 sin 2π L x .The advantage of this example is that all calculations can be performed explicitly (see appendix).
To perform the further discussion we return to equation (3) and consider the force f (x) acting on the particle at its random position as an additional random force f(t) = f[x(t)], which however may be correlated with the noise.Important is, that this force is oscillating in space, and at long times these oscillations take place at a wave vector 2π/L which is much larger than the inverse characteristic scale for the considerable change of the envelope function, being of the order of 1/ √ D * t α .The consideration of such external force corresponds to the substitution This corresponds to the additional contribution 2 v 0 MH(t) ´t 0 f(t ′ )H(t − t ′ )dt ′ in the second and to three additional contributions in the third term on the r.h.s.Now we look at the orders of magnitudes of the corresponding terms, comparing them to the order of magnitude of R 2 (t) ∼ t α .Since H(t) ∼ t α−1 , the first term is subleading: H 2 ∼ t 2α−2 , and 2α − 2 < α for all 0 < α < 2. The contribution of f to the second term is subleading as well, provided f is bounded, |f| ⩽ f 0 : The order of magnitude of this contribution before averaging is For H(t) ∼ t α−1 this contribution behaves as t 2α−1 , and 2α Therefore for subdiffusion this term is subleading as well, so that at long times only the last term survives, so that The correlation function ⟨κ(t ′ )κ(t ′ ′ )⟩ contains four contributions: ⟩, the last one the same as in the force free case and decaying as (t ′ − t ′ ′ ) −α .The corresponding integral of this contribution then gives the term growing as t α , as in the force-free case.Fast equilibration within the single well leads to a fast decorrelation of the values of f(t) = f[x(t)] at different times; the easiest estimate can be obtained for the ⟨f(t ′ )f(t ′ ′ )⟩ and ⟨ζ(t ′ )f(t ′ ′ )⟩ terms with t ′ ′ > t ′ .The behavior under time reversal may be guessed from the standard time-reversal argument in equilibrium (since f (t) equilibrates).The time-dependences of all correlation functions of the type ] corresponds then to a convergent sum of decaying stretched exponentials stemming from single Fourier-components of the Boltzmann factor, where the series is convergent also for t = 0 (note that the only difference of the case of cosinusoidal potential from any other bounded periodic potential is that the coefficients b n can be found explicitly).The contribution of the nth harmonic is a rapidly (as a stretched exponential) decaying function of time, and the long-term behavior of any correlation function of the type considered is dominated by the lowest, first harmonics.Therefore all additional correlation functions containing f (x) are in our case integrable, and their Laplace-transform tends to a constant for s → 0.
Using the properties of stationarity and symmetry under exchange of the time arguments in our, overall equilibrium, noise, we can assume that ⟨κ(t ′ )κ(t ′ ′ )⟩ = F(|t ′ − t ′ ′ |) and then come to the form which is analogous to equation ( 6): The last integral differs from the one in equation ( 7) in that we have in the Laplace domain Noting that, as we have shown, for s small F(s) ≃ K(s) we see that the overall behavior at long times will be the same as without the force.Note that our argumentation is based on the fast local equilibration argument, which is not proved and for which one only has a numerical evidence, but it puts this fast equilibration assumption in connection with the insensitivity of the subdiffusion coefficient to the external force.We note that our discussion can be generalized to potentials built from several Fourier-components (although explicit calculations are unfeasible in this case), and to higher dimensions as well.
Having looked at this simple approximate representation for subdiffusion, one can immediately figure out, under what conditions the external potential still can matter: • A periodic random potential can still strongly influence the motion if it is unbounded, or formally diverges (V 0 → ∞) like in different types of arrangements of solid obstacles, • It can still matter if the external potential (quasi-periodic or random) contains strong lowfrequency components (e.g. in long-range correlated random potentials), • It can also matter for memory kernels which are integrable in time, i.e. in diffusive and superdiffusive cases.
For α = 1 our considerations lead to F(s) ̸ = K(s) and one can await the asymptotic behavior which differs from the one without the force by a prefactor, and use the corresponding equations for obtaining a self-consistent approximation for D * akin to the effective-medium theory (here the renormalization of H(t) is necessary to match the long-time behavior).For superdiffusion, the existence of the external potential, even if the fast local equilibration hypothesis holds, will lead to the change of the power-law in the long-time behavior; presumably it should destroy superdiffusion transforming it into a diffusive motion.
Having this in mind, we numerically simulate one of these situations (particle's subdiffusive motion in a 2D array of solid obstacles, relevant for applications discussed in the Introduction).The discussion of other cases is beyond the scope of the present work.Our main finding is that the long-time dependence of the MSD in the array of impenetrable of scales not influenced by the obstacles' size, just like it was not influenced by the magnitude of a bounded periodical potential in a one-dimensional periodic situation.This is the main finding of the present work.

Subdiffusion in an array of solid obstacles
In what follows we report on the results of simulation of subdiffusion in a periodic 2D array of solid obstacles of approximately circular shape.

A method of numerical simulation
For simulating the subdiffusion we use the Markovian embedding as proposed in [30,37] by solving a system of Langevin equations with n numbering the Laplace modes, α = 1, 2 corresponding to Cartesian coordinates and ξ n,α being the Gaussian white noises with ⟨ζ n,α (t)⟩ = 0 and ⟨ξ n,α (t)ξ m,α (t Here we consider 2D case, which can be clearly visualized; the generalization on 3D case is straightforward and should not introduce new features due to the structure of the system (17).This corresponds to a bath of N Ornstein-Uhlenbeck processes (overdamped bath oscillators) coupled to the particle.As discussed in [30,37], the memory kernel of the corresponding GLE is  (12) this gives 4/π for the diffusivity along each independent coordinate projection, and D * = 8/π for the total simulated 2D system.We chose ν 0 = 10 3 and found C ≈ 1.2874.The initial conditions correspond to equilibrium distribution of the initial coordinates and velocities, namely We integrate the system of equation ( 17) by a forward Euler-Maruyama method.The small time steps required by the Euler integration are not a disadvantage, because they allow us to use a trick for describing the particles' reflections on hitting the obstacles.
Equation ( 17), especially the equation for v, being the Newtonian equation of motion of a tracer particle in the force field corresponding to a sum of the external force and forces due to elastic coupling to the bath oscillators, does literally hold for potentials whose gradients (forces) do not diverge.In our model, this is not the case: while motion between the obstacles corresponds to U = 0 and is force-free, the force on hitting the obstacle diverges, preventing direct integration.The mechanical situation is however extremely simple since hitting an obstacle the particle experiences a specular reflection.This means that we have to know the normal to the obstacle's contour at a point of reflection and to transform the particle's velocity in such a way that its normal component changes the sign and its tangential component stays unchanged.The trick we use to calculate the position of the reflection point, and the coordinate and velocity of the particle after the collision is as follows.
The array of approximately circular obstacles is generated in the following way: we start from generating the 'background field' where S is an empirical correcting summand for the case of a finite n max , and consider the boundaries of the obstacles to correspond to level lines of this function.We then define which function is zero outside of the obstacle and −1 inside of it.The value of n max is chosen to have, for given s, the obstacles which are 'visually round.'Since a large number of summed terms leads to significant slowing down of calculation, we have a trade-off between an approximately circular shape of the obstacles and the computation speed.In our simulations we use n max = 5 and S = 0.55.It is worth noting that for such a small number of harmonics there is no simple dependence between a 'reasonably circular' obstacle's shape and s.We used the following mean radii of the obstacles: R = 0.25, 0.3, 0.4, 0.45 measured in units of the closest distance between centers of obstacles (the respective lines are oriented diagonally) which correspond to s = 0.95, 0.87, 0.42, 0.3.Using larger radii requires larger numbers of components.The packing fraction (portion of the surface unavailable for the motion) is approximately 0.2 for the smallest radius of the obstacles and approximately 0.64 for the largest one (0.81 of the maximal packing fraction for a square lattice under which the obstacles start to overlap).Therefore, the cases we discuss are far enough from the ones when the narrow escape problem through a narrow strait could be of importance.This problem however may also be considered by a dedicated choice of the background field.The outer normal to the boundary of the obstacle is given by the normalized gradient A(x h ) = −∇F(x, y)/|∇F(x, y)| of the field F at the hitting point x h = (x, y).The integration step ∆t is chosen such that flying through the obstacle during a single integration step is virtually impossible.Each main step of integration is subdivided into ten substeps τ = ∆t/10.Staring from the beginning of the ith step at t i we calculate approximate positions x k = x(t i ) + v(t i )kτ , k = 1, . . ., 10, and the corresponding values of f k = f(x k ).The values of f k are zero if the particle is outside of the obstacle and −1 if the trajectory crosses to inside of it.The time from the beginning of the step to hitting the boundary is then ∆t ′ = ∆t + τ 10 k=1 f k and the position at which the particle hits the obstacle is Now we can update the position, where V is the velocity after reflection.The normal component of this velocity is and has the opposite sign to the normal component before the reflection.The tangential component stays unchanged by the reflection.The normal component before the reflection had a direction opposite to those of A and therefore the velocity after the reflection is Now we can update the particle's velocity Finally, the bath variables are updated with Ξ being a Gaussian random variable with zero mean and unit variance and σ = √ 2ν n η n k B T.

Results
Figure 2 shows two examples of single trajectories of the particle in an array of approximately circular obstacles with the average radii R = 0.25 and R = 0.45, which bound the set of examples considered further.The separation between the centers of obstacles is equal to L = 1 in the diagonal direction, the particle starts from the point (0, 0) with equilibrated initial conditions for the velocity, and the maximal time of the motion is t max = 10.This time interval covers the initial ballistic regime and a part of the transient to the long-time subdiffusion regime, where interactions of the random walker with obstacles are of importance, cf figure 3, where the much longer trajectories, up to t max = 10 4 were used for plotting the MSD.
The two panels correspond to different modes of motion: In the case of loose packing one can see relatively long periods of free motion and rare events of encountering with obstacles.On the contrary, in the case of closely packed obstacles, the path represents a tangle of parts of a trajectory concentrated in-between impenetrable areas.Nevertheless, as figure 3 shows, the asymptotic behavior at longer times is the same in both cases.The MSDs calculated based on ensembles of 10 4 realizations for arrays of obstacles of different radii, together with the MSD for the unconstrained subdiffusion, are shown in figure 3 on double logarithmic scales.All curves show a universal ballistic behavior at short times, nonuniversal intermediate-time behavior sensitive to the radius of the obstacles, and the long-time behavior which is again universal, and the same as in the absence of the obstacles.
In the absence of the obstacles, the overall behavior corresponds to a simple and relatively fast crossover from the initial ballistic to the final subdiffusive behavior.In the presence of the obstacles, this crossover is broadened and develops into a transient domain which for largest obstacles' radius stretches over three orders of magnitude in time.The behavior in the transient domain is non-universal but follows a regular pattern: the larger the radius of the obstacles, the earlier the initial ballistic behavior is abandoned, the lower the MSD at a given time is, and the slower is the approach to a final subdiffusive asymptotics.The cases we discuss are still relatively far from the narrow escape regime, which sets on when the obstacles' radii approach R = 0.5 from below.A hint onto a possible additional universal behavior of the transient in the narrow escape case is given by our simulation results for the largest packing fraction considered, corresponding to R = 0.45.Here, the MSD in the intermediate-time domain shows a practically diffusive behavior with R 2 (t) ∼ t over two orders of magnitude in time.On one hand, this is not a wonder in view of the trajectories' shape as illustrated in the right panel of figure 2, where the trajectory spends most of the time captured between the scatterers, and leaves the void between them in a random direction.On the other hand, it is remarkable that multiple scattering events do not lead to the full decorrelation of the motion, temporal memory wins, and the long-time asymptotics of the MSD does not change.Small 'steps' or 'jumps' in the MSD curve at times corresponding the beginning and the end of the diffusive transient highlighted by the dashed line in figure 3 may hint onto additional transient regimes in the narrow escape limit, but these need for a dedicated investigation.The narrow escape variant of the problem however may also be considered within our approach by a dedicated choice of the background field and considerably longer simulation times.

Conclusions
Our work is devoted to a discussion of the subdiffusive particle motion (as described by generalized Langevin equation with a power-law memory kernel) in a periodic array of impenetrable obstacles in two dimensions.
In one dimension, it is known that a bounded periodic potential does not change the asymptotic behavior of subdiffusion under homogenization: the long-time behavior of the MSD in the presence of potential is the same as in the absence of the potential; the influence of the potential is bounded to an intermediate time domain.A qualitative analysis shows that this property may break down if the potential is unbounded or long-range correlated.In the present work we simulate subdiffusion in a periodic 2D array of impenetrable obstacles, corresponding to an unbounded potential, for the case when the obstacles do not overlap and the infinite motion between them is possible.We show that also in this case the long-time behavior is universal, independent on the packing density of the obstacles, and the same as when the obstacles are absent, i.e. that under homogenization, the obstacles do not matter.Our computational method of solving the generalized Langevin equation based on a combination of Markovian embedding as proposed by Goychuk and modeling impenetrable obstacles with 'a background field' trick allows for easy implementation of specular reflection from curved surfaces and allows for circumventing a usually computationally complicated task of modeling the motion in nonsmooth, unbounbded external potentials.Our simulations were restricted to the case when the straights between the obstacles were relatively large (i.e.outside of the narrow escape variant of the problem), but the method can be generalized also to such situations as well as to other situations with impenetrable boundaries of complex topology, like the ones appearing in modeling anomalous diffusion in cells and tissues.
The revealed universality shows that one has to be extremely cautions when addressing the subdiffusion detected from the MSD dependencies obtained in macroscopic experiments with brain tissues, e.g. by studying the spread of concentrated drops of fluorescent markers.In this case, in contrast to the case of Brownian motion (normal diffusion), it is hard to draw definitive conclusions on the microscopic topology of the brain's ECS (e.g. its porosity, tortuosity, etc) from macroscopic data.

Appendix. The case of cosinusoidal potential
As we already stated, for the special case of a cosinusoidal potential U(x) = − f0L 2π cos 2π L x , corresponding to the sinusoidal force f(x) = −f 0 sin 2π L x , calculations can be performed explicitly.
Putting down with a 1 = π + arcsin f f0 and a 2 = − arcsin f f0 , and substituting the expression for p(x, t) given by equation ( 13), we get one easily recognizes an integral sum for a Gaussian integral, so that the overall sum converges to 2. Therefore, at long times the distribution of the force is with C being the normalization constant, and is equal to the equilibrium distribution of the force in a single well of the potential.One can also consider the two-point distribution.
Let us now consider the time-dependences of correlation functions of the type with Y(t) being ζ(t) or f[x(t)], with t ′ > t, and estimate the averaging over the realizations of the trajectories using the PDF of displacements while leaving averaging over the initial conditions explicit (and for later).Using the quasi-equilibrium distribution, equation ( 13) and putting formally t = 0 with x(t = 0) = x 0 we get: where the average is taken over the initial distribution of x 0 which might be reduced to the first elementary cell of the potential, and Z = ´L 0 e − U(x 0 ) k B T dx 0 .Now we denote the Boltzmann factor by g(x) and note that for a periodic potential this is a periodic function of x and can be expressed by a Fourier series.For the cosinusoidal potential the series reads a sinus-Fourier series with decaying non-negative coefficients (note that for real z ⩾ 0 one has I n+1 (z) < I n (z), see e.g.[39]).Returning to our integrals in x for the correlation function we see that this is represented by the sum of Fourier harmonics of the Gaussian envelope function.This one is easy to obtain representing sine as a difference of complex exponentials and noting that the characteristic function of a Gaussian centered at x 0 is given by exp(i x 0 k − σ 2 k 2 /2) with k now being a wave vector: Therefore the contribution of the nth harmonic is a rapidly (as a stretched exponential) decaying function of time, and the long-term behavior of any correlation function of the type considered is dominated by the lowest, first harmonics.Therefore all additional correlation functions containing f (x) are in our case integrable, and their Laplace-transform tends to a constant for s → 0. Since the approach described above addresses the fast decay of individual high harmonic components of the characteristic function, it can be straightforwardly generalized to smooth bounded potentials possessing several periodic components, i.e. represented by the Fourier series.Moreover, one can also 2D regular soft potentials given by superposition of two periodic functions each depending on one coordinate only.

ORCID iDs
Eugene B Postnikov  https://orcid.org/0000-0001-7904-1881Igor M Sokolov  https://orcid.org/0000-0002-4688-9162 with the Ms term being a leading one.The dependence of the correlation function of the noise (being the property of the bath only) on the mass of the test particle makes the superballistic situation non-generic, if possible at all.

Figure 1 .
Figure 1.(A) Cosinusoidal potential stated by equation (14) with parameters L = √ 2, f 0 = 1.(B) The probability distribution function of the random walk with α = 0.5 in such potential for t = 10 4 (solid blue curve) superposed with the regression line of this curve by a parabola (dashed red curve) on a logarithmic scale.
and the noise correlation function fulfills equation (2).The reasonable choice of ν n is ν n = ν 0 b −n , where ν 0 and b are the high-frequency cutoff and the dilation parameter, respectively, and η n = C(α, b)ν α n .The constant C(α, b) has to be found empirically to provide the best approximation for K(t) ≃ K 0 t −α for given N over the longest time span.We use the decade scaling b = 10 with N = 16 components.Our simulation parameters are: α = 0.5 (relatively deep in the subdiffusion domain), M = 0.1, k B T = 1 and K 0 = 1.According to equation

Figure 2 .
Figure 2. Examples of trajectories of the particle in arrays of solid obstacles with radii R = 0.25 (A) and R = 0.45 (B) measured in units of length equal to the shortest distance between centers of the neighboring obstacles in the diagonal direction.

Figure 3 .
Figure 3.The particles' MSD in an arrays of solid obstacles of different radii as indicated in the legend (the unit length is the shortest distance between centers of the nearest obstacles).The exponent of anomalous diffusion is α = 0.5, see text for details.The green dotted straight line represents the asymptotic long-time behavior for the case when the scatterers are absent; the segment of a dashed gray straight line indicates that under high packing fraction the behavior in the intermediate time domain is approximately diffusive.