Cusp of Non-Gaussian Density of Particles for a Diffusing Diffusivity Model

We study a two state “jumping diffusivity” model for a Brownian process alternating between two different diffusion constants, D+>D−, with random waiting times in both states whose distribution is rather general. In the limit of long measurement times, Gaussian behavior with an effective diffusion coefficient is recovered. We show that, for equilibrium initial conditions and when the limit of the diffusion coefficient D−⟶0 is taken, the short time behavior leads to a cusp, namely a non-analytical behavior, in the distribution of the displacements P(x,t) for x⟶0. Visually this cusp, or tent-like shape, resembles similar behavior found in many experiments of diffusing particles in disordered environments, such as glassy systems and intracellular media. This general result depends only on the existence of finite mean values of the waiting times at the different states of the model. Gaussian statistics in the long time limit is achieved due to ergodicity and convergence of the distribution of the temporal occupation fraction in state D+ to a δ-function. The short time behavior of the same quantity converges to a uniform distribution, which leads to the non-analyticity in P(x,t). We demonstrate how super-statistical framework is a zeroth order short time expansion of P(x,t), in the number of transitions, that does not yield the cusp like shape. The latter, considered as the key feature of experiments in the field, is found with the first correction in perturbation theory.


Introduction
The emergence of non-Gaussian features for the positional probability density function (PDF) of particle spreading, denoted P(x, t), in a disordered environment is a common attribute that arises in many different physical and biological systems. Specifically, a tent like shape of the PDF, in the semi-log scale, together with a linear time dependence of the mean square displacement (MSD) appear for diffusion in glassy system [1], biological cells [2][3][4][5][6][7], and colloidal suspensions [8][9][10][11]. This tent shape, sometimes fitted with a Laplace distribution P(x, t) ∼ exp(−C|x|) with C a constant, suggests that the decay of the PDF is exponential. This feature is becoming a more frequent observation for the spreading of molecules. Phenomenological approaches are diffusing diffusivity models, in which non-Gaussianity is obtained by coupled stochastic differential equations with random diffusion coefficients [7,[12][13][14][15][16][17][18][19][20], and path integrals formalism for Brownian motion in the presence of a sink [21]. More recently, theoretical frameworks describing this behavior emerged from continuous time random walk (CTRW) approaches employing large deviations theory [22][23][24] and microscopical models like molecular dynamics of tracer particles in polymer networks [25,26] and interacting particles with fluctuating sizes [27][28][29], the so-called Hitchhiker model [28].
While in some of the systems the non-Gaussian behavior disappears when the measurement time is made long enough, the short time tent-like decay of the PDF seems to be a universal phenomenon [22]. It is then natural to ask if there is some sort of universality that can be deduced for the temporal limit of short times. Within the diffusing diffusivity models for the large x limit, exponentially decaying propagators have been observed by employing a dichotomous process for the diffusivity [13]. The latter model consists of a "fast" and a "slow" phases, each one with a diffusion coefficient D + and D − , respectively [13,30]. Furthermore, the appearance of a cusp at small displacements also has been reported in different diffusive approaches like the Sinai model [31], employing the quenched trap model [32][33][34][35][36][37] or spatial dependence in the diffusivity [38], within the Lévy-Lorentz gas model [39] and using the fractional Fokker-Planck equation [40]. It is important to notice that the cusp found in [31][32][33][34]36,[38][39][40] is within the context of anomalous diffusion in the MSD sense, and those presented in [35][36][37] are for normal diffusive systems.
It is worth mentioning that several systems in nature exhibit, or can be reduced to, a dichotomous process. Examples of two state systems include nuclear magnetic imaging to measure the diffusion of heterogeneous molecules [41], diffusion in glassy materials [1], blinking quantum dots [42,43], diffusion in single molecules tracking experiments [4,7], and protein conformational dynamics [44]. Other approaches for analyzing two state systems were also devised over the years; see heterogeneous molecular transport [41], telegraphic noise [43], Lèvy Flights [45], and CTRW models [1,22].
In this work, we deal with a two state jumping diffusivity model with equilibrium initial conditions, i.e., we assume that the process started long before the measurement began. The long measurement time behavior of the positional PDF for this model is Gaussian and is independent of the specifics of the waiting times at the different diffusive states. A rather unexpected result is achieved for the opposite temporal regime. We obtain that the behavior in the limit of the short measurement times and the shape of the positional PDF of the molecule spreading in the two state jumping diffusivity model attains a cusp or a general tent-like shape. Our result is based on the statistics of the temporal occupation fraction of the diffusivity states, the latter is defined as the time spent in state D + over the total measurement time. The Gaussian behavior in the long measurement time is dictated by the δ-function shape of the distribution of this temporal occupation fraction, a feature that is solely based on the ergodic properties of the system. We show that, in the limit of short measurement time, the distribution of the temporal occupation fraction attains a uniform distribution that leads to the mentioned cusp behavior of P(x, t). The uniformity of the occupation fraction is a general result in the sense that it does not depend on the statistics of the waiting times in the two states, and the latter can be arbitrary. The non-Gaussian behavior of P(x, t) for short measurement times is similarly general as the Gaussian behavior of the propagator for long times. We then show that our approach reproduces the results of a specific representative system with exponentially distributed waiting times.
Our manuscript is organized as follows: in Section 2.1, we introduce the jumping diffusivity model and the initial conditions utilized in this work. In Section 2.2, we develop our theory for the statistics of the occupation time in the short measurement time limit, for which the PDFs of the waiting times in states D ± are rather general. The obtained behavior of the occupation fraction is used in order to describe the non-Gaussian features of P(x, t), i.e., its cusp shape, that is observed in this model. In Section 2.3, we corroborate our previous results for a system with exponentially distributed waiting times. In Section 3.1, we discuss briefly how these theoretical results differ from those found within the super-statistical approach [12,14,21] and further how our approach may be applicable in experiments. Finally, in Section 4, we present a summary of our results, and we discuss briefly recent work of Postnikov et al. [37] who considered a model with quenched disorder, emphasizing the importance of equilibrium initial conditions. The main derivations are given in the corresponding Appendixes.

The Model
We consider a two state renewal model, with a stochastic diffusion field D(t) for a particle in a random medium. The position of the particle is following a diffusion process given by dx(t)/dt = 2D(t)ξ. D(t) ∈ {D + , D − } is a dichotomous model, considering the case when D − < D + and ξ is a standard white noise, i.e., with mean zero, variance one, and delta correlated. As an example of the dynamics of the model, at a given time, the particle follows a pure diffusion process with a diffusion coefficient D + > 0 during a period τ. After this time period has elapsed, the diffusion coefficient jumps and, during the next time interval, the particle diffuses with diffusion coefficient D − . The waiting times at each state D ± are distributed according to a general PDF ψ ± (τ), with mean waiting times τ ± . The subscript ± denotes whether the waiting times are defined for the D + or D − states. In the following, we present the two-state model with D − = 0, while the case with D + > D − > 0 is analyzed in Appendix A. In Figure 1, we show representative trajectories for the position at time t, x(t), while, in Figure 2, we present the same for D(t) and we show the notation we use. t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 8 τ 1 τ 2 τ 3 τ 4 τ 5 τ 6 τ 7 τ 8 Figure 1. Typical trajectory of x(t) given by Equation (2) with D + = 10 (blue regions), and D − = 0 (red regions). For this trajectory, exponential waiting times with τ + = 1 and τ − = 5 were used.
We define T ± as the occupation time in state "±", namely the total amount of time that the process diffuses with D + or D − during t. Jumps between states D + and D − occur at random times t 1 , t 2 , etc., until a final measurement time t and clearly t = T + + T − . The intervals of time between each jump are defined by Figure 2. Then, the occupation times in each state, when started from D + , are explicitly provided by where N is the random number of transitions that were performed between the two states during the measurement time t and k is an integer. The measurement time t and N satisfy t ≥ t N , with t N = τ 1 + τ 2 + . . . + τ N , i.e., the exact time when the Nth jump was performed. The backward recurrence time τ * is defined by τ * = t − t N [46]. Each waiting time τ i follows τ i = t i − t i−1 with i ∈ (1, N). For this particular initial condition, odd values of i in τ i refer to waiting times at D + and even values of i to waiting times during which the diffusion coefficient is D − (see Figure 2). Expressions similar to Equation (1) are also obtained when the process starts from D − , see Equation (A65) . Since the particle is diffusing with a constant diffusion constant D + for time τ 1 , when starting from D + , the position x(τ 1 ) is simply x(τ 1 ) = √ 2D + τ 1 ξ 1 , where ξ 1 is a zero mean Gaussian random variable with ξ 2 1 = 1. When at the state with diffusion constant D − = 0, the particle is not moving, therefore where all ξ i are independent zero mean Gaussian random variables that satisfy ξ 2 i = 1. By using Equation (1) and exploiting the properties of summation of independent Gaussian variables, we obtain that the position at general time t is provided by when D + > 0 and D − = 0. Equation (2) holds irrespective of the state at t = 0. We see that the particles' position is a product of two independent random variables, the square root of the time staying at the state D + times a standard Gaussian random variable. Alternating process for the diffusivity, starting from the state '+' and N = 2k + 1. For the case of equilibrium initial conditions exposed in Section 2.2, for N = 1, τ 1 works as the forward recurrence time with PDF Equation (10).
In the following, we consider a situation in which the process has started long before the measurement began, i.e., at t = 0, the process was already running for a very long time. In this way, the measurement begins from an initial condition in which the system is in equilibrium, meaning that the probability to start from D + is τ + /[ τ + + τ − ] and accordingly the probability to start from [46,47]). For this set-up, the PDF of the occupation time T + , f t (T + ), is determined by the contribution to start from D + and the contribution to start from D − , yielding where f ± t (T + ) is the PDF of T + for measurement time t, given that the process has started from ±. Since D − = 0, Equation (2) dictates that the positional PDF, provided that the system has occupied the state with D + for a time T + , is given by The propagator of the system is obtained via integrating over all possible values of the occupation time T + , whose PDF is f t (T + ) Equation (3), yielding Likewise, we can work with the temporal occupation fraction, which is defined by p + = T + /t with 0 ≤ p + ≤ 1. In this case, the positional PDF for a specific value of p + follows and the propagator is obtained similarly to Equation (5), but using the PDF of p + , which we denote by g t (p + ), Since the properties of P(x|T + ) or P(x|p + ) are known, the task of computing the propagator completely depends on our ability to calculate the PDF of T + or p + . In the following section, we address this problem.

The General Case: Arbitrary Distribution of Waiting Times
Two regimes of the process are of special interest. The long and the short limits of the measurement time t. The two different limits involve different considerations when computing the PDFs of the occupation time (T + ) and fraction (p + ). We first handle the regime of small t and then we treat the t → ∞ limit.

Short Time Regime
The PDF of the occupation time T + is defined by Equation (3). We condition on the number of transitions N, and each term f ± t (T + ) is provided by where Q ± t (N) is the probability to perform exactly N transitions during t when the process started at ±. f ± t (T + |N) is the PDF of T + when exactly N transitions were performed (during t), and the process has started from ±. This conditional probability is obtained by counting the number of trajectories of temporal span t that started from the ± state and performed exactly N transitions, out of the total number of trajectories that started from the ± state and for which the diffusion spent a total time T + at this state. Utilizing Equation (8), we rewrite Equation (3) as Since we consider a renewal process, the expression for Q ± t (N) is known in the Laplace space [42], asQ ± s (N) = L{Q ± t (N)} = ∞ 0 Q ± t (N) exp(−ts) dt, for any generalψ ± (s) = L{ψ ± (τ)}. Concretely, Q ± t (N) is obtained by taking into account all the possibilities to perform N jumps up to time t N < t, and no additional jumps during the backward recurrence time τ * . This sums up to a convolution of N + 1 random variables. It is important to notice that, since we assume equilibrium initial conditions, τ 1 , which is measured from t = 0, is only a part of a full renewal event and is termed the forward recurrence time. The PDF of τ 1 for the ± state, f ± eq (τ 1 ), is provided by (see [46]) and in the Laplace space L{ f ± eq (τ 1 )} = (1 −ψ ± (s)) τ ± s. This initial condition stems from the equilibrium of the underlying process, in which we do not have a jump at the initial time (t 0 = 0 in Figure 2). In the literature [13,42,46,48,49], the case where the renewal process starts at t = 0 is called ordinary or non-equilibrium, and as we will see below, by following our approach, this does not yield any universal features for P(x, t), hence the assumption of an equilibrium process is important in our methodology, (see discussion about non-equilibrium initial conditions in Appendix B).
The probability of not performing any jumps during τ * is equivalent to the probability of obtaining a waiting time τ N+1 > τ * , i.e., 1 − τ * 0 ψ ± (τ) dτ. Eventually, by implementing the initial equilibrium condition, we obtain In all the equations above on the right-hand side, we have a multiplication of functions in the Laplace space, this implies convolutions as we transform from s to t. The first term in the multiplication on the right-hand side of Equation (11) obviously stems from the equilibrium initial condition under study. We assume that the PDF of the waiting times is analytic for τ → 0, thus we can express ψ ± (τ) as [22,23] with A ± ≥ 0 an integer number. As an example, consider the case with exponential waiting times, i.e., ψ ± (τ) = ψ(τ) = exp(−τ/ τ )/ τ , namely the waiting times at the D ± states are identically distributed. Its analytic expansion is ψ The analyticity of ψ ± (τ) Equation (12) is a very mild demand that covers a wide range of sojourn times distributions. Since we are interested in the small t limit, the corresponding behavior in the Laplace space is found for s → ∞, where the leading terms ofψ ± (s) are [22] ψ ± (s) ∼ For the mentioned example with exponential waiting times,ψ(s) ∼ 1/[ τ s]. Using Equation (13) forQ ± s (N), we obtain that, in the s → ∞ limit, corresponding to the short time limit, which is at the focus of our interest We see that the leading terms for allQ ± s (N) with N > 1 are of the order 1/s γ with γ > 2. Thus, in the small t limit, terms with N > 1 contain contributions that scale like t γ−1 and are negligible with respect to the N ∈ {0, 1} cases. Therefore, only the first two Q ± t (N)s are taken into account, i.e., This is an expected result, as, for short times, only contributions from a single transition and zero transitions are important. By calculating Q ± t (N), we advanced towards obtaining the behavior of the PDF of T + , according to Equation (9), in order to complete this mission, one needs to compute the relevant contributions of f ± t (T + |N) in the t → 0 limit. First, we see that the conditional distribution f ± t (T + |0) depends only on the starting state. There are only two types of trajectories that have performed 0 transitions, i.e., for all the time, they have been either at D + or at D − . Consequently, The calculation of f ± t (T + |N) is obtained by conditioning over the first event. If starting from the + state, the process will spend a time τ 1 at this state before jumping to the − state. τ 1 can attain any value 0 ≤ τ 1 ≤ T + and for the remaining time t − τ 1 the process has to perform one transition less. In general, without regarding the initial conditions of the problem, an integration over all possible τ 1 's provides the relation with N + 1 = N, and B + a normalization factor. For instance, for N = 1, we have that t 0 ψ + (τ 1 )/B + dτ 1 = 1, which stems from the fact that we consider only trajectories of time span t. The corresponding formula for Since we are assuming equilibrium initial conditions, the ψ ± in the N + 1 element of the iterative forms (Equations (19) and (20)) must be replaced by f ± eq (Equations (10)). As was already noted above, only the N = 0 and N = 1 are of interest in the small t limit, then, according to Equations (17), (20), and (10), we get for N = 1 Using the small time approximation of ψ ± (τ) Equation (12), in Equations (21) and (22), we obtain that, independently of the starting state, The 1/t dependence comes from the integral factors in Equations (21) and (22), all the other terms in the numerator and denominator simply cancel out. See Appendix B for a complementary derivation of Equation (23) using the definition of the joint PDF of T + and N. Gathering Equations (15)- (18), and (23) in Equation (9), we find that The PDF of the occupation fraction is obtained trivially from Equation (24) by changing variables to The third term in Equations (24) and (25) is uniform, i.e., terms which are independent of T + or p + , and this is the first main result of this paper. All the additional terms and contributions to the PDF of p + only introduce terms that depend on higher orders of t and are thus negligible in the small t limit. This means that, for equilibrium initial conditions, regardless of the exact form of ψ ± (τ), the PDF of p + (Equation (25)) is always uniform for 0 < p + < 1. This general uniform behavior of the PDF of the occupation fraction is applicable for an extremely large class of waiting times PDFs ψ ± (τ). As a remark, the connection between the conditional PDF of T + , f ± t (T + |N), and the joint PDF of T + and N, f ± t (T + , N) is discussed in Appendix B.1. In the following, it is shown that this uniformity leads to universal features of the propagator in the limit of small t. In Sections 2.3.1 and 2.3.2, we treat particular examples (with exponential waiting times) that are exactly tractable, without any simplifications or assumptions. The results agree perfectly with the general form in Equation (25). It is important to notice that our approximations affect only the form of g t (p + ) and do not affect P(x|p + ). This allows us to obtain the behavior of P(x, t) for any −∞ < x < ∞, as is shown below.

P(x, t) for Arbitrary Waiting Times
In order to obtain the positional PDF for small t, we combine Equations (6), (7), and (25), which, after integration, gives Considering After substituting in Equation (26), it turns into with We can see that, in Equation (27), there is a linear dependence on |x| in the vicinity of x = 0. This means that for short enough measurement times the PDF of x will always have a tent like shape, irrespective of the distributions ψ ± that were chosen (see Figure 3 below). Only the mean sojourn times affect the shape. This is a general result for the short time regime, and it is based on the general fact that the PDF of the temporal occupation fraction is uniform for 0 < p + < 1. Concretely, at short times when |x| is small, the decay of P(x, t) will always resemble an exponential one. For large |x|, the form of P(x, t) must be Gaussian, due to the fact that this limit is determined by the instances when no transition to D − was ever made and the transport is controlled by diffusion with D + . However, if we only look at the particles that have moved, i.e., we get rid of the delta function at x = 0 in Equation (26). We can relate these dynamics with some experiments which condition the measurements on the movement of the particles. This procedure is called population splitting; see [50,51]. Technically, if D − > 0, the cusp is not found; however, as long as D − /D + << 1, the tent like shape will be found; for further details, see Appendix A.

Long Time Regime
In the limit t −→ ∞, the PDF of the temporal occupation fraction g t (p + ) follows a different but also a general form. As mentioned, we are focusing on the case where both ψ ± have finite first moments, τ ± > 0. In the long time limit, ergodicity is satisfied, namely the equivalence of ensemble and temporal averages are attained. Particularly, in this case, the ensemble average of the occupation fraction at D + is equal to the temporal average which is defined by the fraction of average waiting times at D + and D − , i.e., p + = τ + /[ τ + + τ − ] (see Appendix F). Thus, in the long time limit, g t (p + ) converges to a δ-function, Since ergodicity prevails, by using Equation (28) in Equation (7), the positional PDF gets the form In the long time limit, the positional PDF given by Equation (29) represents a Gaussian propagator with an effective diffusion coefficient Since τ + /[ τ + + τ − ] < 1, and the effective diffusion coefficient is always smaller compared with D + . Indeed, this slow-down is an expected result due to the portion of the time that the particle spends in the state with D − = 0 and basically not moving during this period.
Both cases fit with Equation (26) (red and blue solid lines) with a tent like shape. In both normalized histograms at x = 0, there is a peak representing the Dirac delta function in Equation (26). Right: P(x, t) for t = 30 and waiting times uniformly distributed (green triangles) with the same parameters as above and for gamma distributed waiting times (orange squares) with We employed the last set of parameters in the gamma distributed waiting times in order to avoid an overlapping between curves. P(x, t) converges to the Gaussian statistics Equation (29) (green and orange solid lines). In all the presented cases, D + = 10 and D − = 0 were used.

Simulations
The two general limits of g t (p + ) Equations (25) and (28) produce two different prevailing distributions of P(x, t) Equations (26) and (29). In Figure 3, we compare analytical formulas Equations (26) and (29) (solid lines) with simulations of two different state modelsone with uniform distributed waiting times τ ∼ U(0, 5) for D + and τ ∼ U(0, 10) for D − and with t = 1 (red triangles) and t = 30 (green triangles), such that τ + = 2.5 < τ − = 5. Here, the notation τ ∼ U(a, b) means that τ has a uniform distribution with a and b the minimum and maximum values, respectively. In addition, the other with gamma distributed waiting times, such that τ ∼ Gamma(k, θ). The latter notation denotes that τ has a gamma distribution with k its shape parameter and θ the corresponding scale parameter. In this case, the PDF follows particularly the PDF of the gamma distribution Equation (30) implies a cumulative distribution function F(τ) = γ(k, τ/θ)/Γ(k), with γ(x, y) the incomplete gamma function and Γ(x) the standard gamma function. For the latter case, we used τ ∼ Gamma(0.5, 8) at D + and τ ∼ Gamma(0.5, 12) at D − , for t = 2 (blue squares) and t = 30 (orange squares), with τ + = 4 < τ − = 6. As we can see in the short time regime, for uniform and gamma distributed waiting times (red triangles and blue squares), P(x, t) has a tent shape for short displacements, and it agrees with Equation (26), joined with a peak at x = 0 due to the Dirac delta function in Equation (26). For t = 30 (green triangles and orange squares), each case of P(x, t) converges to Gaussian statistics (Equation (29)).
The cusp we have found for small |x| implies that we may approximate the distribution on a small scale with a Laplace like distribution, P(x, t) ∼ exp(−C|x|). However, clearly this does not hold globally for large x, see Figure A2 in Appendix C. Still within the interval of short displacements, due to the presence of the delta peak at the origin, we expect for this span a considerable contribution on the normalization of P(x, t). Particularly, we find that the area underneath the curve for the case of uniformly distributed waiting times (red line) in Figure 3 in the left panel has a value of 0.88 for x ∈ (−4, 4). In addition, the corresponding area within the same figure, but, for gamma distributed waiting times, the (blue curve) has a value of 0.89 for x ∈ (−8, 8).

Exponentially Distributed Waiting Times
In this section, we obtain g t (p + ) for a specific distribution of waiting times, but using different methods, which let us corroborate the validity of our general approach described above. We analyze the case of exponential waiting times in states with D + and D − , each waiting time following a PDF given by We show first the case of a two state system with the same mean waiting times and then investigate the complimentary case. In Appendix D, we analyze both cases for non-equilibrium initial conditions, e.g., a system starting from D + .

Equal Mean Waiting
Times τ + = τ − Let us consider a system with τ + = τ − = τ . We know that the temporal fraction occupation p + and T + can be related to the difference of occupation times defined by [46]. In this section, we analyze the double The Laplace transform of ψ(τ) in Equation (31) is given by L ψ(τ) =ψ(s) = 1 1+ τ s . Substitutingψ(s) in Equation (32), we obtain In Appendix E, an analytical expression for the PDF of S t is found, i.e., inverse Laplace transform of Equation (33) is performed (see Equation (A44)). Then, remember that the temporal occupation fraction in the plus state p + is related to the difference of occupation times as S t = 2p + t − t. We can employ Equation (A44) for obtaining the PDF of p + , which is given by A similar expression for a system with non-equilibrium initial conditions (always starting from D + ) is found in Appendix D. Expanding Equation (34) in the short time limit t −→ 0, i.e., t << τ , Equation (34) can be approximated by a uniform distribution For τ + = τ − = τ , Equation (35) agrees with Equation (25) obtained by the general approach of Secton 2.2. In the left panel of Figure 4, we show the short time approximation of g t (p + ) (Equation (35)) compared with the general formula in Equation (34); it is evident that both results agree perfectly. In the right panel of Figure 4, we show Equation (34) for short and long measurement times. g t (p + ) evolves from a uniform distribution to a peaked distribution centered at its mean value p + = 1/2 (see Appendix E for a deduction of the central moments of g t (p + )).

Positional Distribution Function
An analytical expression for the positional distribution function P(x, t) (given by Equation (7)), with g t (p + ) provided by Equation (34), can be deduced by using the series representation of the modified Bessel functions, I ν (y) = ∞ ∑ k=0 ( y 2 ) 2k+ν /[k!Γ(ν + k + 1)]. The integration in Equation (7) yields with 1 F 1 (a; b; z) the confluent hypergeometric function of the first kind. Nonetheless, in a short time limit, we can use the uniform approximation of g t (p + ) (Equation (35)), and then Equation (7) provides which agrees with the results obtained above in Equation (26), when τ + = τ − and for t −→ 0, since, in that limit, exp(−t/ τ ) ∼ 1 − t/ τ . Particularly for x = 0 and taking x −→ 0, Equation (37) yields a tent shaped propagator described by with K 2 = (5t − τ )/[16 τ √ π(D + t) 3 2 ] and in concordance with Equation (27). On the other hand, within this short time limit, for large displacements x −→ ∞, the two terms between curly braces in Equation (37) cancel each other, and only the first term in Equation (37) is left (when x = 0). This is due to the expansion of 1 Then, Equation (37) can be approximated by This Gaussian behavior of P(x, t) at the tails is expected. The large |x| limit is dominated by trajectories for which no transitions to D − were performed and a pure diffusion process with D + occurs. When t >> τ , ergodicity is satisfied and, therefore, the system on average visits the two states the same amount of time. Namely, the ensemble average of p + is equal to the corresponding fraction of the average waiting times. In this case, when τ + = τ − , the occupation fraction is concentrated at p + = 1/2. Thus, the PDF of p + is represented by the delta function Substituting Equation (40) in Equation (7), we recover Gaussian statistics for the displacements In Figure 5, we present the two different limit distributions for P(x, t) in the short time limit t = 0.1 (red circles) and t = 0.5 (blue crosses) Equation (37) and the Gaussian limit for t = 5 (orange circles) and t = 10 (green crosses) Equation (41), for the normalized variable z = x/ √ t. As we can see, the displacements for short times follow a tent shape (black solid line) and a Gaussian one in the long time limit (magenta solid line).

Different Mean Waiting Times τ + = τ −
Relaxing the assumption of equal mean waiting times for exponentially distributed sojourn times in the model, we have that τ + = τ − , with waiting times following Equation (31). As mentioned, for equilibrium initial conditions, the PDF of T + is given by Equation (3) Then, the different terms of the PDF of T + in Equation (3) are provided in Laplace space, by [42,47,49] Summing up Equations (42) and (43) according to Equation (3), we obtain, for exponentially distributed waiting times, Taking the double inverse Laplace transform of Equation (44) with respect to u ⇔ T + and s ⇔ t and changing variables to p + = T + /t, we obtain the PDF for p + (see details in Appendix F) For the case when τ + = τ − = τ , Equation (45) recovers Equation (34) obtained by the methods reported in [52,53]. The case of non-equilibrium initial conditions is shown in Appendix D.
In the short time regime, strictly speaking when t << τ ± , by expanding Equation (45) for t −→ 0, g t (p + ) can be approximated by the uniform distribution As mentioned above, Equation (25) that was deduced for general PDFs of waiting times, encloses the particular case of Equation (46). For the uniform approximation of g t (p + ) (Equation (46)), the positional PDF (Equation (7)) is which agrees with the general case described by Equation (26). Similar to Section 2.2, in the limit t −→ ∞, the PDF of the occupation fraction g t (p + ) follows Equation (28). In addition, the PDF of the displacements in the long time regime is given by Equations (7) and (29), recovering Gaussianity.
In Figure 6, we show g t (p + ) for exponential waiting times with τ + = 1 and τ − = 5, in the left panel, we compare the uniform approximation of Equation (46) (black asterisks) with the full solution Equation (45) (red solid line), observing an excellent agreement. In the right panel of Figure 6, the behavior of g t (p + ) (as provided by Equation (45)) is displayed. As we can see, it starts with a uniform distribution for short times and then it evolves to a peaked distribution centered at p + = τ + /( τ + + τ − ) = 1/6. As shown in Appendix D, for non-equilibrium initial condition, the PDF of p + is still uniform within the short time regime. See also Appendix B for other similar cases. Finally, in Figure 7, we show the corresponding positional spreading for the normalized variable z = x/ √ t. As we can see in the short time, t = 0.1 (red circles) and t = 0.5 (blue crosses) P(z, t) (given by Equation (47)) attain a tent-like shape. In the long run, t = 20 (orange circles) and t = 30 (green squares) P(z, t) have a Gaussian distribution given by Equation (29). 7. For a system with, τ + = 1 and τ − = 5, P(z, t) in semi-log scale, with z = x/ √ t. For short times t = 0.1 (red circles) and t = 0.5 (blue crosses), P(z, t) is represented by Equation (47) (black solid line) with a tent like shape. For large times t = 20 (orange circles) and t = 30 (green diamonds), P(z, t) converges to the Gaussian statistics Equation (29) (magenta solid line). In all the cases, D + = 10 and D − = 0 were used. Compared with Figure 5, in this case, the Gaussian curve is above the tent curve, contrary to the case with equal mean waiting times. This is because the coefficient of the Gaussian curve Equation (29) is bigger compared with the weight of the delta peak in Equation (47). In Figure 5, we have the opposite, and the weight of the corresponding delta function in Equation (37) is bigger compared with the Gaussian Equation (41).

Super-Statistics
We have found that at x = 0, P(x, t) exhibits a cusp. A mathematically similar non-analytical behavior is found using an approach called super-statistics [12,14,21,54], which was used to explain laboratory observations. This framework postulates that the distribution of diffusion constants in the system is exponential, namely P(D) = exp(−D/ D )/ D for D > 0 and D the average diffusivity. Then, the diffusion follows a Gaussian process with a random D. This approach gives Here, on the right-hand side, we have the Laplace PDF, which was used by Laplace in 1774 [55] to describe his linear law of errors [56]. In addition, as in our case, within the super-statistics method, we see in Equation (48) a non-analytical behavior since P(x, t) ∼ C 1 − C 2 |x|, for small x and with C 1 , C 2 constants. Our work does not support the Laplace law, see Equations (26) and (27). However, maybe more importantly, the whole approach presented in this manuscript differs from the super-statistical approach in the following way. In our model, we have two diffusion constants, D + and D − = 0 (see Appendix A for the case when D − = 0). Hence, the PDF of diffusion constants is P(D) = aδ(D) + bδ(D − D + ), with a, b ≥ 0. It follows that the super-statistical approach predicts that the diffusing packet P(x, t) is a sum of a delta function corresponding to non-moving particles and a Gaussian packet describing the movers. Thus, when the non-moving particles are excluded, we have perfect Gaussian behavior. This is actually correct, to leading order, for very short times. Thus, the super-statistical approach gives the correct t −→ 0 behavior but fails to predict the main issue (in our opinion), and that is the cusp on x = 0. To explore the non-analytical behavior, one needs to go to the next order terms in the expansion to include paths with a transition between states. Then, as we have shown, the equilibrium initial condition yields a uniform distribution of the occupation fraction Equation (25). It is this fact that brings the non-analytical behavior in the final result for P(x, t) Equation (27), graphically represented by a "tent" see Figures 3,5 and 7. It follows that the exponential conspiracy in which distribution of diffusion constants is exponential is not a necessary condition for a cusp like behavior of P(x, t). We further remark that the non-analytical behavior is found also in the context of normal diffusion in [35][36][37] and within the anomalous one at [31][32][33][34]36,[38][39][40].

Time Average MSD
We note that, in single molecule experiments, the time average mean squared displacement (TAMSD) is used in many cases to estimate the distribution of diffusion constants [2,5,6]. Since time averages are recorded over a finite measurement time, the time average fluctuates. Hence, we have naturally a distribution of the estimator for the diffusion parameters. In addition, the aforementioned two delta peak distribution of D, i.e., on D + and on D − , is expected to be smeared out. This topic was extensively studied in a wide variety of models [57,58].
We now investigate the fluctuations of the time averaged diffusivities in a two state model and their implications in the distribution of diffusion coefficients obtained from real experimental data. For a further analysis of the time average diffusivity within a two state system, see [59,60].
We note that, in different single particle tracking experiments with non-Gaussian propagators, the recorded distribution of the diffusion coefficient D (obtained by means of TAMSD analysis) is relatively broad and peaked close to the origin [2,5,6]. Those experimental distributions of D are typically fitted by exponential [6] or gamma [2] distributions. Within the two state model, the diffusivity takes only two possible values D − or D + , but the respective TAMSD analysis gives values of D around D − and D + [59]. The average D is given by D = (D + τ + + D − τ − )/( τ + + τ − ). Thus, how different is the distribution of the diffusivities, extracted via TAMSD techniques, in a two state model compared with the one present in single molecule experiments? As we show next, this will be determined by the values of D ± and τ ± . In Figure 8, we show the distribution of the diffusion coefficients obtained by means of TAMSD analysis for D + = 10, D − = 0 and exponentially distributed waiting times. We show two different cases, the first one with the same mean waiting times τ + = τ − = 1 (see red boxes). In addition, the second one with different mean waiting times, such that τ + = 1 and τ − = 5 (see blue boxes). From the linear plots of the TAMSD versus the lag time estimates of D were extracted. We show two cases, the first for a system with the same mean waiting times τ + = τ − = τ = 1 (red boxes). In addition, the PDF of D for a system with different mean waiting times with τ + = 1 and τ − = 5 is also shown (blue boxes). For the system with the same mean waiting times, the average diffusivity found in the simulations is D = 4.98, and, for the case of different mean waiting times, we have D = 1.69. In both cases, we used t = 1000 and 1000 trajectories.
As we can see in Figure 8, when the difference between the diffusion coefficients is large, as in our case D + = 10 > D − = 0, P(D) is relatively broad. Nonetheless, for the case with τ + = 1 and τ − = 5, the peak of P(D) is closer to the origin compared to the case with τ = 1.
This difference between mean waiting times in each state is the second factor that determines the shape of P(D). For instance, when this difference is such that τ + < τ − , it is straightforward that the more the process spends in the state "−", the more the observed values of D will be closer to D − . In this latter case, the distribution of D is peaked close to the origin since D − < D + . Thus, we can say that, when the differences between the diffusivities (and the mean waiting times) in the different states are pronounced, i.e., D − << D + and τ + << τ − , P(D) in the two sate model resembles the distributions found in single molecule experiments [2,5,6].

Conclusions
From symmetry of the density of spreading particles P(x, t) = P(−x, t), we expect an analytical expansion of the propagator as P(x, t) ∼ K 1 − K 2 x 2 + . . ., with K 1 , K 2 constants. Instead, in the two state model handled throughout this work, we get an expansion that is linear in |x|, see Equation (27). This is a non-analytical expansion graphically represented by a tent like structure, see Figures 3, 5 and 7. As mentioned above, Laplace in 1774 considered a similar non-analytical PDF, P(x) = exp(−|x|)/2 for −∞ < x < ∞ [55,56]. However, the expression we find is clearly non-exponential, see Equation (26). Furthermore, for large x, we get a Gaussian behavior for P(x, t). It should be noted that a non-analytical behavior is found only if D − = 0, see Appendix A for further details. In practice, we may approach the non-analytical features of P(x, t), as D − is getting small.
Recently, a very general theory was developed for the non-Gaussian spreading of packets of particles. Using a CTRW framework, it was shown that, for any analytical PDF of waiting times, for large x limit P(x, t) ∼ exp(−C|x| ln |x|), with C a constant [22]. In the former model, we thus find exponential tails for large x, while, here, the anomaly, i.e., the cusp or tent like feature of P(x, t), comes from the small x limit.
Recently, Postnikov et al. [37] investigated a model of diffusion in a quenched disordered setting, where the diffusive field is spatially varying. They showed that equilibrium initial conditions play a major role stating: "within the class of models with quenched disorder, the Itô model under equilibrium conditions is the only promising candidate for the description of Brownian Non Gaussian diffusion (BnG)." Note that here the definition of BnG means a model or system where the MSD is increasing linearly for all times and the propagator is non-Gaussian. Our model uses a time dependent diffusivity, and we showed that equilibrium initial conditions are indeed a key requirement. Here, we note that BnG does not imply a cusp, and vice versa. Namely, we may find a system where the MSD is increasing linearly in time, for the entire span of time, with or without a cusp for P(x, t) at x = 0. The main focus of our work is the presence of a cusp for P(x, t). Regarding the behavior of the MSD, it can be shown that, when equilibrium initial conditions are applied, T + = ( τ + t)/[ τ + + τ − ], for all times t (see Appendix G). Then, by Equations (A62) and (A68), the MSD is provided by for any time t. Thus, if the process starts from equilibrium, the MSD grows linearly for all times and we have BnG. Nevertheless, we would like to emphasize that our model is exhibiting BnG, but specifically P(x, t) has a cusp only if D − = 0, and practically when To summarize, we emphasize that we have shown, by means of the statistics of the temporal occupation, that there is a universality for the PDF of the temporal occupation fraction in a two state model. For PDFs of waiting times with finite first moments, g t (p + ) can be approximated by a uniform distribution following Equation (25). This leads to tent like decaying propagators (Equation (26)) similar to those found in many experimental systems. We corroborate our results by solving analytically a two state system with exponentially distributed waiting times. We have shown that, either for short or long times, the distribution of displacements P(x, t) has a general form, either a "tent" or a Gaussian bell curve. These two endpoints of the positional PDF are independent of the actual form of the distribution of waiting times. The crucial point within our framework is the generality of the behavior of the PDF of the occupation fraction p + , being a uniform distribution for short times and a delta peak for long times. The former was found for a system with equilibrium initial conditions. We note that, for certain types of non-equilibrium initial conditions, we can still get a uniform PDF for the fraction occupation time; however, this is not generic (see details in Appendix B). Therefore, the non-Gaussian features are readily present in our model within the short time regime, and regardless of the specifics of the waiting times.
Mathematically, we presented an expansion in terms of the number of transitions from state + to − and backwards. Naturally, for very short times, the leading contribution to the packet comes from the paths with zero transitions, and then the packet is simply a sum of two Gaussian curves with diffusion coefficients D + and D − . However, we showed that, by going to next order terms in the expansion, namely considering the paths with a single jump, we get the cusp like shape, found in the limit D − → 0. Thus, the whole effect is achieved by using a perturbation approach obtaining the leading order correction to the trivial behavior. Put differently, a widely popular super statistical approach is found to miss one of the main issues of the field, namely the cusp in P(x, t). A super-statistical approach [14,54] uses a distribution of diffusivities, which in our model is a sum of two delta functions, at D − = 0 and D + . This does not give the cusp, as it is merely the zeroth order of the perturbation theory developed here.  Data Availability Statement: All the data sets obtained by numerical simulations or data analysis are available from the corresponding authors upon request.

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

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. A Two State Model with D + > D − > 0
When D + > D − > 0, the process for the displacements becomes with ξ 1 and ξ 2 each i.i.d. Gaussian variables. In this case, the form of the conditioned PDF is given by Then, the marginal distribution for the displacements follows g t (p + )dp + . (A3) Appendix A.1. P(x, t) for Arbitrary Waiting Times As we did in Section 2.2.1.1, using the general forms obtained above, i.e., Equations (9), (17), (18), (23), and (28), we can analyze P(x, t) in the short and long time limits. Figure A1. Distribution of displacements P(x, t) obtained by simulations of a two state system with D + > D − > 0 and gamma distributed waiting times τ ∼ Gamma(3, 1) at D + and τ ∼ Gamma(6, 1) at D − following Equation (30). We compare with Equation (A4) (solid lines) with t = 0.5, τ + = 3, Substituting Equation (25) in Equation (A3), we get In Figure A1, we compare Equation (A4) (in solid lines) and P(x, t) obtained by simulations of a two state model for a fixed value of D + and different values of D − , such that D + > D − . In all cases, we used gamma distributed waiting times τ ∼ Gamma(3, 1) for the state with D + and τ ∼ Gamma(6, 1) for the state with D − (the gamma distribution is defined by Equation (30)). As we can see when D + >> D − , e.g., D + = 10 and D − 0.1 (red triangles), the PDF of the displacements at small values of x has a non-Gaussian peak; thereafter, for large values of x, it follows a Gaussian distribution. When the values of D − approach D + (cyan squares and magenta circles), P(x, t) is fully described by Gaussian statistics even in the short time limit.
Appendix A.1.2. Long Time Regime In the long time limit, the PDF of temporal occupation fraction is provided by Equation (28), then, according to Equation (A3), the PDF of the displacements is determined by Thus, the Gaussian limit is also restored.
Appendix A.2. P(x, t) for Exponentially Distributed Waiting Times with τ + = τ − In the short time regime, we can use the uniform approximation Equation (46) in Equation (A3), and the distribution for the displacements yields For the long time regime, we use g t (p + ) provided by Equation (28); then, according to Equation (A3), the PDF of the displacements follows Gaussian statistics described by Equation (A5).
with the distribution of jumps defined by [46] with 1 (a,b) (t) the indicator function, such that it is equal to one if t ∈ (a, b) and zero if t / ∈ (a, b). The average · is over all the the values of τ i 's, with i ∈ {1, 2, . . . , N + 1}.In addition, f ± t (T + , N) is the joint probability of T + and N, which satisfies [46] with t N = τ 1 + . . . + τ N and the average · defined as above.
Let us find Equations (A8) and (A9) and therefore Equation (A7) for the case N = 1. It is important to notice that, for the case of equilibrium initial conditions such as the one handled in Section 2.2, when N = 1, the corresponding average on τ 1 is given by the forward recurrence distribution Equation (10). Following Equation (A8), and taking the Laplace transform defined asQ ± which is already the result shown in Equation (11).
For the joint distribution f ± t (T + , 1), following Equation (A9) and taking the double N)dtdT + , after performing the corresponding integrals in the case we started from "+", we get Following the same procedure for the case when the process started from "−", we obtain Next, we show the connection of the joint PDf Equation (A9) with the uniform distribution of the occupation times Equation (24). Letf eq s (u, 1) be the double Laplace transform of f eq t (T + , 1), i.e., the joint PDF of T + and one single jump, starting from equilibrium. Clearly, the former followŝ f eq s (u, 1) Using Equations (A11) and (A12) in Equation (A13), we get One of the key features of our paper is found when we consider both u and s to be large. This corresponds to the short time limit, when T + and t are of the same order. Then, we may useψ + (s + u),ψ − (s) −→ 0 in Equation (A14), yielding tô Equation (A15) is easy to invert, and we find in the short time limit In order to analyze the joint, we have to deal with the analytical expression of ψ ± (τ) defined by Equation (12). For instance for N = 2, after substituting Equation (13) in Equations (A23) and (A25), we have that the dominant term isf + , respectively. Inverting the double Laplace transform, in each case, it gives a positive power of T + and therefore, with respect to t, e.g., A ∓ ≥ 0 is a positive integer number, this correction term for t −→ 0 (T+ −→ 0) is negligible, and also the remaining terms in f ± t (T + , N) with N > 2. We conclude that, for the case of equilibrium initial conditions, the uniformity in the short time limit, for the PDF of the occupation/fraction time, is always preserved, as long ψ ± (τ) is analytical.

Appendix B.1. Non-Equilibrium Initial Conditions
Still, by employing the joint distribution f ± t (T + , N), it can be shown, as follows that, for non-equilibrium initial conditions, when ψ ± (τ) in the Laplace space is approximated byψ ± (s) ∼ 1/s for short times, it can lead to a uniform distribution in the occupation time and therefore to a tent shape in P(x, t).
In the case of non-equilibrium initial conditions, either starting just from D + or from D − , the averages over τ 1 , within Q ± t (N) Equation (A8) and the joint distribution f ± t (T + , N) Equation (A9), are no longer given by f ± eq (τ 1 ) Equation (10). Now, in the non-equilibrium case, these corresponding averages are performed using the waiting time PDF ψ ± (τ 1 ). Following the same procedure as above, for the case N = 1, the double Laplace transform of the joint distribution f ± t (T + , 1) yields tô For a system with non-equilibrium initial conditions, in order to recover the uniform distribution in the PDF of T + , it is enough to ask that, for large s (short times), the PDF of the waiting times in the Laplace space followŝ Substituting Equation (A27) in Equation (A26), we have thatf ± s (u, 1) ∼ 1/[(s + u)s]. This implies in the real space, for 0 < T + < t, that the joint distribution follows f ± t (T + , 1) ∼ 1, and therefore, because of Equation (3), we also have a uniform distribution for T + . As an example of a model in whichψ ± (s) goes as Equation (A27) and for non-equilibrium initial conditions ψ ± (τ) = f ± eq (τ 1 ) as expected, we have the case in which the PDF of waiting times is the sum of two exponential functions, e.g., ψ ± (τ) For this latter case, sinceψ ± (s) satisfies Equation (A27), following the same analysis as above, we find that the joint distribution f ± t (T + , 1) is uniform. In Appendix D, we show that, for a system with non-equilibrium conditions and exponentially distributed waiting times with equal and different mean values, the distribution of the occupation/fraction time is also uniform. This is mainly because, for exponentially distributed waiting times, the forward recurrence time distribution in equilibrium f ± eq (τ 1 ) Equation (10) is equal to ψ ± (τ 1 ), as in the non-equilibrium case. Furthermore, for exponentially distributed sojourn times, Equation (A27) is also satisfied, and its PDF in the Laplace space for s −→ ∞ followsψ ± (s) ∼ 1/[ τ ± s]. By using Equation (3), this gives the uniform distribution shown in Equations (A35) and (A40).
Forψ ± (s) given in Equation (A27), the correction terms when N ≥ 2, following the same analysis as in the case of equilibrium initial conditions, yield to elements of the form , which are negligible for t −→ 0 ⇐⇒ T + −→ 0. Thus, they do not contribute in the PDF of the fraction occupation time.

Appendix C. P(x, t) from Simulations with Uniform and Gamma Distributed Waiting Times within the Complete Range of x
Following Figure 3 in the left panel in which, for short time and displacements, the cusp of P(x, t) is displayed. Now, from simulations (with the same parameters as above) of a two state model, with uniform (red triangles) and gamma (blue squares) distributed waiting times. In Figure A2, we show P(x, t) in semi-log scale but for the whole span of x.
In each case, we compare the normalized histogram of the simulation data with the short time analytical formula Equation (26), finding a perfect agreement. As we can see, the cusp is located at the origin, and, for large displacements, Gaussianity is recovered. Figure A2. Distribution of displacements P(x, t) in semi-log scale, obtained from simulations, of a two state system with uniform and gamma distributed waiting times within the short time limit and displaying the whole span of x. P(x, t) for uniformly distributed waiting times is shown in red triangles. In addition, the case of gamma distributed waiting times is shown in blue squares. We employed the same set of parameters as those used in Figure 3 in the left panel. Both cases fit with Equation (26) (red and blue solid lines). state with D + is 1 and the probability of starting from the state with D − is 0. The PDF of T + then satisfies With f + t (T + , N) given by Equation (A9) explicitly for this case, we have [46] f + t (T + , 2k Thus, using Equation (A30) for summing over all the values of N in Equation (A28), we get For exponentially distributed waiting timesψ(s) = 1/(1 + τ s), by substitutingψ(s) in Equation (A31), we get that the double Laplace transform of Equation (A28) is given bŷ By the same procedures used in Appendix F, the inversion of the double Laplace transform of Equation (A32) yields Employing the identity I ν (y) = (y/2) ν 0F1 (; ν + 1; y 2 /4) [61] and changing variables, we obtain the PDF of the occupation fraction, which follows g t (p + ) = δ(1 − p + )e − t τ + t τ I 0 2t τ p + (1 − p + ) By taking the series expansion of Equation (A34) in the limit t −→ 0, the PDF of p + can be approximated by Therefore, for 1 > p + > 0, the PDF of p + follows a uniform distribution (see the left panel of Figure A3), as in the case of equilibrium initial conditions (Equation (35)). Figure A3. Left: g t (p + ) Equation (A34) for τ = 1 and t ∈ {0.1, 0.5, 1, 2, 5, 10} and non-equilibrium initial conditions (starting from state "+"). The uniform approximation of g t (p + ) Equation (A35) for t = 0.1 is shown in black circles. Right: g t (p + ) Equation (A39) for τ + = 1, τ − = 5 and t ∈ {0.1, 0.5, 2, 5, 10, 20} and non-equilibrium initial conditions (starting from state "+"). The uniform approximation of g t (p + ) Equation (A40) for t = 0.1 is shown in black circles.
The inverse Laplace transform of Equation (A38) is obtained by the same procedures explained above and, in Appendix F, eventually g t (p + ) = δ(1 − p + )e − t τ + + t τ + I 0 2t p + (1 − p + ) We recover Equation (A34) when τ + = τ − = τ . In the short time limit, Equation (A39) follows as well a uniform distribution for 1 > p + > 0, see the right panel of Figure A3. In this case, the PDF of p + is given by Thus, for exponentially distributed waiting times, for equilibrium and non-equilibrium conditions, in the short time regime, the PDF of the occupation fraction is always uniform. This feature is only valid for exponentially distributed waiting times, and it is not necessarily fulfilled for other distributions of waiting times. P(x, t) for the specific case of non-equilibrium initial conditions and exponentially distributed waiting times is the double Laplace transform (t ⇔ s and T + ⇔ u) of the PDF of T + Equation (3) In the second line of Equation (A62), we have employed the linearity of · , and the properties of independent standard normal random variables, i.e., ξ 2 i = 1 and ξ i ξ j = 0 (with i, j ∈ {1, 2} and i = j). Thus, now we just have to find T + . In order to do that, we start from the definition of average occupation time Withf t (u) = ∞ 0 f t (T + )e −uT + , and T + ⇔ u Laplace conjugates. Now, for equilibrium initial conditions, the PDF of the occupation time is given by with f ± t (T + , N) the joint PDF of the occupation times at D + and the number of jumps between states during t, once started from D ± . When starting from D + and having N = 2k + 1 or N = 2k jumps, the occupation time in each case satisfies Equation (1). In the case when the initial state is at D − , we have T + = τ 2 + τ 4 + . . . + τ * i f N = 2k + 1, with τ * = t − t N , the backward recurrence time. The definition of the joint PDF f ± t (T + , N) is already given in Equation (A9). In addition, its double Laplace transformf ± s (u, N) = Now, for obtainingf s (u), we compute the double Laplace transform of Equation (A64) and then we sum Equations (A22), (A25), and (A66) for all values of N. Thereafter, we compute the derivative off s (u) with respect to u and its corresponding limit when u −→ 0. Following algebraic simplifications, we yield lim u→0 d duf s (u) = − τ + τ + + τ − 1 s 2 . (A67) For obtaining the average occupation time Equation (A63), we invert Equation (A67) with respect to s, having Finally, substituting Equation (A68) in Equation (A62), we get Equation (49), which indicates that the MSD is linear with respect to t, for any value of time t.