<figure><img src="image_1.png" /></figure>

HAL
open science

# Deriving sub-diffusion equations

Benoît Perthame, Min Tang

► To cite this version:

Benoît Perthame, Min Tang. Deriving sub-diffusion equations. 2025. hal-05185833

HAL Id: hal-05185833

https://cnrs.hal.science/hal-05185833v1

Preprint submitted on 25 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Deriving sub-diffusion equations

Benoît Perthame*

Min Tang†

July 24, 2025

## Abstract

Sub-diffusion equations are used in a large range of applications including fluids, plasma physics and biology. Their mathematical analysis is advanced even if a much larger literature addresses super-diffusions. The goal of this paper is to provide the microscopic mechanism and rigorous derivation of sub-diffusions when the waiting time distribution of particles follows an age-structured equation and jumps occur at each renewal. The major difficulty to recover sub-diffusions, unlike normal diffusions, is that the assumption of long waiting time implies lack of integrability for the age equilibrium. This prevents to establish strong a priori estimates. Here, the Laplace transform plays the role that Fourier transform plays for the more traditional case of fast diffusions.

2020 Mathematics Subject Classification. 35R09, 35B40, 44A10, Keywords and phrases. Structured population dynamics; Multiscale analysis; Caputo derivatives;

## 1 Introduction

Sub-diffusion, characterized by a mean squared displacement (MSD) scaling as $\langle x^2(t) \rangle \propto t^\alpha$ with $0 < \alpha < 1$, has been observed across diverse physical and biological systems [19, 26, 25, 16]. In living cells, sub-diffusive transport is prevalent: lipid granules in fission yeast exhibit anomalous diffusion in the cytoplasm [25], insulin granules in MIN6 cells show heavy-tailed waiting time distributions during microtubule binding/unbinding [21], and hydration water on protein surfaces undergoes sub-diffusion due to trapping in heterogeneous potential landscapes [22, 23]. For a comprehensive review of sub-diffusion in natural systems, see [18].

Mathematically, sub-diffusion is often modeled using fractional calculus, introducing nonlocal time derivatives to account for memory effects. The Caputo Riemann-Liouville fractional derivative, defined as

$$ {}_0^{RL}D_t^\alpha f(t) = \partial_t \int_0^t \frac{f(\tau)}{(t-\tau)^\alpha} d\tau, $$

naturally arises in sub-diffusion equations such as

$$ {}_0^{RL}D_t^\alpha u - \operatorname{div}(\kappa(u)\nabla u) = f, $$

<footnote>*Sorbonne Université, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris.
Email: Benoit.Perthame@sorbonne-universite.fr</footnote>

<footnote>†School of Mathematical Sciences and Institute of Natural Sciences, MOE-LSC, CMA-Shanghai, Shanghai Jiao Tong University, China. Email: tangmin@sjtu.edu.cn</footnote>
where $\alpha$ governs the degree of sub-diffusion [7]. This formalism captures memory through the power-law kernel $(t - \tau)^{-\alpha}$.

A prevalent microscopic model for sub-diffusion is the continuous-time random walk (CTRW), where particles alternate between jumps and trapping events with stochastic dwell times. When the trapping time distribution follows a power law decay $\psi(t) \sim t^{-(1+\alpha)}$ ($0 < \alpha < 1$) for $t \gg 1$, the ensemble-averaged MSD exhibits sublinear scaling $\langle x^2(t) \rangle \propto t^\alpha$. This non-ergodic behavior arises because longer observations uncover deeper traps, leading to a breakdown of ergodicity [19].

Prior works have linked CTRW models to sub-diffusion equations. For instance, [3] showed equivalence between their Green's functions, while [2] derived the fractional Fokker-Planck equation from a generalized master equation for CTRW, connecting power-law trapping kernels to fractional time derivatives. In CTRW models, the power-law trapping distribution is typically assumed, which implies that the trapped particles are always in equilibrium. However, a fundamental question remains: *What microscopic mechanisms generate power-law trapping distributions?*

In this paper, we address this gap by proposing an **age-structured model** for particle trapping dynamics. Age-structured models, widely used in demography, ecology, epidemiology [24, 10, 8], track populations partitioned by age $a$ (here, trap residence time). We model the probability density $u(a, t)$ of particles trapped for time $a$ at time $t$ using the age structured (McKendrick-Von Foerster) equation [15]:

$$ \frac{\partial u(a, t)}{\partial t} + \frac{\partial u(a, t)}{\partial a} = -d(a)u(a, t), \qquad (1) $$

where $d(a)$ is the age-dependent escape rate. When $d(a) = \frac{\alpha}{1+a}$ ($0 < \alpha < 1$), the escape probability decreases with trap residence time, leading to a power-law waiting time distribution $\Psi(a) \sim a^{\alpha-1}$ for large values of $a$.

By coupling this age-structured trapping model with spatial jumps, we rigorously derive the sub-diffusion equation and identify scaling conditions under which the microscopic and macroscopic descriptions coincide. Notably, we establish a correspondence between the initial conditions of the age-structured model (e.g., initial trap residence times) and those of the fractional sub-diffusion equation, addressing a critical gap in the literature. This connection provides mechanistic insights into sub-diffusive behavior in complex systems. An age-structured kinetic model was used in [11] to derive a space-time fractional diffusion equation, but the derivation is formal, and the initial condition was not identified.

On the one hand, the mathematical literature on problems involving **Caputo derivatives** and partial differential equations is extensive, spanning topics from the porous medium equation [1] to mathematical biology and medicine [4, 14] (see also the references therein). Sobolev regularity in $L^p$ spaces has also been explored (cf. [9]), quasi-linear problems have been investigated [27], and Besov regularity has been established. However, these works do not address the derivation of such equations from finer-scale mechanisms. On the other hand, equations combining age-dependent escape rates with kinetic equations have been derived from particle scattering studies [5, 17]. Various limiting models have also been investigated: for instance, the authors of [12] and [13] analyze the normal diffusion limit and super-diffusion limit, respectively. A consistent assumption in their work is that the equilibrium in Eq. (1)—specifically, $u_{eq} = e^{\int_0^a d(a')da'}$—is integrable. In our analysis, a major difficulty is that the sub-diffusion regime emerges precisely when this equilibrium is non-integrable. This corresponds to long waiting times between jumps and also hinders the establishment of strong a priori estimates. Finally, to our knowledge, this is the first rigorous derivation of a sub-diffusion equation from a microscopic model via asymptotic analysis, with appropriate initial conditions explicitly identified. This regime is examined in [6] under a different scaling, yielding a Hamilton-Jacobi equation, while self-similar
solutions are constructed in [4]. Our approach differs from these works, and—analogous to the use of Fourier transforms for super-diffusion—we employ Laplace transforms to derive the sub-diffusion equation.

The organization of this paper is as follows: Section 2 introduces the notations and demonstrates how to apply the Laplace transform to prove the decay rate for large time. Section 3 establishes the power-law decay rate for an age-structured model in the large-time regime, while Section 4 derives the sub-diffusion equation when spatial jumps are modeled by a nonlocal kernel. The corresponding scales of the age-structured jump model are identified, and the contribution from the initial condition is shown to emerge as a source term in the sub-diffusion equation.

## 2 Laplace transform and time decay

As a preparation to the case with space, we introduce notations and material concerning the Laplace transform and its use for proving large time decay rates. We consider functions $u(t)$, $t \in \mathbb{R}$ satisfying $u(t) = 0$ for $t \le 0$ and $u(t) \ge 0$ for $t > 0$. The Laplace transform of $u(t)$ is denoted by $\hat{u}(s)$,

$$ \hat{u}(s) = \mathcal{L}_s u = \int_0^\infty e^{-st} u(t) dt, \quad s \ge 0. $$

Notice that, thanks to our convention for $t \le 0$, $s \hat{u}(s) = \int_0^\infty e^{-st} \frac{d}{dt} u(t) dt$, which can be rewritten as

$$ s\mathcal{L}_s u = \mathcal{L}_s \partial_t u. \qquad (2) $$

We also remark that, for $u \ge 0$,

$$ u \in L^1(\mathbb{R}) \iff \hat{u} \in L^\infty(0, \infty). $$

Moreover,

$$ u(t) \approx t^{\alpha-1}, \quad t \gg 1, \quad 0 < \alpha < 1 \quad \text{implies} \quad \hat{u}(s) \approx s^{-\alpha} \to \infty \quad s \approx 0, $$

because, with the change of variable $t \mapsto \sigma = st$,

$$ \int_0^\infty t^{\alpha-1} e^{-st} dt = s^{-\alpha} \int_0^\infty \sigma^{\alpha-1} e^{-\sigma} d\sigma. \qquad (3) $$

Here since $\alpha \in (0, 1)$, $\sigma^{\alpha-1}e^{-\sigma}$ is integrable on $[0, \infty)$. The reverse implication is not so simple and can only be interpreted in an integral sense as given by Lemma 1.

**Lemma 1** Assume that, $u(t) \ge 0$ for $t > 0$, and that for some $\alpha$ with $0 \le \alpha \le 1$ and constants $K, \bar{K}$, the following holds for $0 < s \le 1$: $Ks^{-\alpha} \le \hat{u}(s) \le \bar{K}s^{-\alpha}$. Then for all $\varepsilon > 0$ and for some constants $\underline{K}_1(\alpha)$ and $\bar{K}_1(\alpha)$ that are independent of $\varepsilon$,

$$ \frac{\underline{K}_1(\alpha)}{\varepsilon} \le \int_0^\infty u(t) t^{-\alpha-\varepsilon} dt, \quad \int_{t \ge 1} u(t) t^{-\alpha-\varepsilon} dt \le \frac{\bar{K}_1(\alpha)}{\varepsilon}. \qquad (4) $$

**Proof.** We prove the second inequality. We first use the immediate control

$$ \int_0^1 s^{-1+\alpha+\varepsilon} \int_0^\infty e^{-st} u(t) dt ds = \int_0^1 s^{-1+\alpha+\varepsilon} \hat{u}(s) ds \le \bar{K} \int_0^1 s^{-1+\varepsilon} ds = \frac{\bar{K}}{\varepsilon}, $$
and thus we conclude, with the change of variable $s \mapsto \sigma = st$, that

$$ \int_{1}^{\infty} u(t) \int_{0}^{1} s^{-1+\alpha+\varepsilon} e^{-st} ds dt = \int_{1}^{\infty} u(t) t^{-\alpha-\varepsilon} \int_{0}^{t} \sigma^{-1+\alpha+\varepsilon} e^{-\sigma} d\sigma dt \leq \frac{\overline{K}}{\varepsilon}. \quad (5) $$

This gives the upper bound because $t \geq 1$ and $\int_{0}^{1} \sigma^{-1+\alpha+\varepsilon} e^{-\sigma} d\sigma \geq \frac{e^{-1}}{1-\alpha-\varepsilon} > \frac{e^{-1}}{1-\alpha}$.
For the lower bound, the same calculations give, on the one hand

$$ \int_{0}^{1} s^{-1+\alpha+\varepsilon} \int_{0}^{\infty} e^{-st} u(t) dt ds \geq \frac{K}{\varepsilon}, $$

and on the other hand, similar equality as in (5) gives

$$ \frac{K}{\varepsilon} \leq \int_{0}^{\infty} u(t) t^{-\alpha-\varepsilon} \int_{0}^{t} \sigma^{-1+\alpha+\varepsilon} e^{-\sigma} d\sigma dt \leq \int_{0}^{\infty} u(t) t^{-\alpha-\varepsilon} dt \int_{0}^{\infty} \sigma^{-1+\alpha+\varepsilon} e^{-\sigma} d\sigma. $$

Then $K_1(\alpha)$ is given by $K_1 = \int_0^\infty \sigma^{-1+\alpha+\varepsilon} e^{-\sigma} d\sigma$.

**Remark 2** The conclusion (4) is a statement compatible with the large time decay of $u(t)$. Namely, when $u(t) \approx t^{\alpha-1}$, (4) is satisfied. Notice also that the same conclusion follows with the weaker assumption for $s$ such that

$$ \frac{K}{\varepsilon} \leq \int_{0}^{1} s^{-1+\alpha+\varepsilon} \hat{u}(s) ds \leq \frac{\overline{K}}{\varepsilon}. $$

### 3 Age structured equation and time decay

To apply the above considerations about the Laplace transform, we consider the age structured equation:

$$ \left\{ \begin{array}{l} \partial_t u(t, a) + \partial_a u(t, a) + d(a)u(t, a) = 0, \quad a \ge 0, t \ge 0, \\ U(t) := u(t, 0) = \int_0^\infty d(a)u(t, a)da, \\ u(0, a) = u^0(a), \end{array} \right. \quad (6) $$

and $u(t, a)$ is extended to $0$ for $t < 0$. Here $d(a) > 0$ and the initial data is supposed to satisfy

$$ u^0(a) \geq 0, \quad u^0(a) \in L^1(0, \infty), \quad M := \int_0^\infty u^0(a) da < \infty. \quad (7) $$

The problem is conservative, i.e.,

$$ \int_{0}^{\infty} u(t, a) da = M, \quad \forall t \geq 0. $$

We define

$$ D(a) = \int_{0}^{a} d(a') da'. \quad (8) $$

When $e^{-D}$ is integrable, it is standard that solutions of (6) converge to a multiple of $e^{-D}$ [15, 20]. When it is not, it remains that $de^{-D}$ is integrable. For this reason, we investigate the equation satisfied by the renewal term $U(t)$.
**Useful formulas.** In order to avoid abstract technicalities, we often use

$$d(a) = \frac{\alpha}{1+a}, \qquad e^{-D(a)} = \frac{1}{(1+a)^{\alpha}}, \qquad \alpha \in (0, 1). \qquad (9)$$

Because $\alpha$ is fixed, we set

$$\Gamma_{\alpha} := \Gamma(1 - \alpha) = \int_{0}^{\infty} \frac{e^{-\sigma}}{\sigma^{\alpha}} d\sigma. \qquad (10)$$

For any given function $P(a) \in L^1[0, \infty)$, the following formula is useful

$$\mathcal{L}_s \partial_t \int_0^t \frac{P(a)}{(t-a)^\alpha} d\tau = \Gamma_\alpha s^\alpha \mathcal{L}_s P. \qquad (11)$$

This formula holds because, integrating by parts in $t$, we can compute

$$\int_0^\infty e^{-st} \partial_t \int_0^t \frac{P(a)}{(t-a)^\alpha} dadt = s \int_0^\infty e^{-st} \int_0^t \frac{P(a)}{(t-a)^\alpha} dadt = s \int_0^\infty P(a) e^{-as} \int_a^\infty \frac{e^{-s(t-a)}}{(t-a)^\alpha} dtda \\ = s \int_0^\infty P(a) e^{-as} \int_a^\infty \frac{e^{-st}}{t^\alpha} dtda = s^\alpha \int_0^\infty P(a) e^{-as} da \int_a^\infty \frac{e^{-\sigma}}{\sigma^\alpha} d\sigma,$$

where the last equality follows from the change of variable $t \mapsto \sigma = st$.

In this paper, we also use a related formula

$$\mathcal{L}_s \partial_t \int_0^t \frac{P(\tau)}{(1+t-\tau)^\alpha} d\tau = s \int_0^\infty \frac{e^{-sa}}{(1+a)^\alpha} da \mathcal{L}_s P. \qquad (12)$$

It is derived exactly as for (11).

**The equation for $U(t)$.** Applying the Laplace transform to the age structured model (6), we find

$$\partial_a \hat{u}(s, a) + (s + d(a)) \hat{u}(s, a) = u^0(a).$$

The solution is given by

$$\hat{u}(s, a) = \hat{u}(s, 0)e^{-D(a)-sa} + e^{-D(a)-sa} \int_0^a u^0(a')e^{D(a')+sa'} da'$$

Thus, the Laplace transform of the renewal term $\hat{U}(s) := \hat{u}(s, 0)$ is determined by

$$\hat{U}(s) = \hat{U}(s) \int_{0}^{\infty} d(a)e^{-D(a)-sa} da + \int_{0}^{\infty} d(a)e^{-D(a)-sa} \int_{0}^{a} u^{0}(a')e^{D(a')+sa'} da' da. \qquad (13)$$

To treat the first term on the right hand side (RHS) of (13), we compute

$$\int_0^\infty d(a) e^{-D(a)-sa} da = -s \int_0^\infty e^{-D(a)-sa} da + \int_0^\infty (s+d(a)) e^{-D(a)-sa} da = -s \int_0^\infty e^{-D(a)-sa} da + 1.$$
Therefore, (13) becomes, using the Fubini theorem and the same argument as before but replacing the lower bound of the integrations from 0 to $a'$,

$$ 
\begin{aligned}
s \int_0^\infty e^{-D(a)-sa} da \hat{U}(s) &= \int_0^\infty u^0(a')e^{D(a')+sa'} \int_{a'}^\infty d(a)e^{-D(a)-sa} dada' \\
&= \int_0^\infty u^0(a) da - s \int_0^\infty u^0(a')e^{D(a')+sa'} \int_{a'}^\infty e^{-D(a)-sa} dada'.
\end{aligned}
$$

When $d(a)$ is given by (9), the above equation can be written as

$$ s \int_0^\infty \frac{e^{-sa}}{(1+a)^\alpha} da \hat{U}(s) = \int_0^\infty u^0(a) da - s \iint_0^\infty u^0(a') \frac{(1+a')^\alpha}{(1+t+a')^\alpha} e^{-st} dt da'. \quad (14) $$

Here the last equality is obtained by the change of variable $t = a - a'$.

Using the formula (2) and (12), we may rewrite (14) as

$$ 
\begin{aligned}
\mathcal{L}_s \partial_t \int_0^t \frac{U(\tau)}{(1+t-\tau)^\alpha} d\tau &= \int_0^\infty u^0(a) da - s \mathcal{L}_s \int_0^\infty u^0(a) \frac{(1+a)^\alpha}{(1+t+a)^\alpha} da \\
&= \int_0^\infty u^0(a) da - \mathcal{L}_s \partial_t \int_0^\infty u^0(a) \frac{(1+a)^\alpha}{(1+t+a)^\alpha} da,
\end{aligned}
$$

and thus

$$ \partial_t \int_0^t \frac{U(\tau)}{(1+t-\tau)^\alpha} d\tau = \delta(t) \int_0^\infty u^0(a) da - \partial_t \int_0^\infty u^0(a) \frac{(1+a)^\alpha}{(1+t+a)^\alpha} da. \quad (15) $$

We finally obtain the equation that $U(t)$ satisfies, which is in the following integral form:

$$ \int_0^t \frac{U(\tau)}{(1+t-\tau)^\alpha} d\tau = \int_0^\infty u^0(a) da - \int_0^\infty u^0(a) \frac{(1+a)^\alpha}{(1+t+a)^\alpha} da. \quad (16) $$

For $t > 0$, differentiating both sides of (16) with respect to $t$ yields (15), whereas both sides of (16) vanish at $t = 0$. The integral on the left-hand side of (16) has no singularity in the integration kernel; hence, it is not a Caputo derivative but still incorporates the memory effect.

From the explicit exact formula (16), we have the time decay property of $U(t)$:

**Proposition 3 (Time decay for age structured eq.)** With assumptions (7) and (9) with $\alpha \in (0, 1)$, the solution of (6) decays as $u(t, 0) = U(t) \approx t^{\alpha-1}$ for $t \gg 1$, in the sense that

$$ \lim_{t \to \infty} \int_0^t \frac{U(\tau)}{(1+t-\tau)^\alpha} d\tau = \int_0^\infty u^0(a) da. \quad (17) $$

In particular, for all $0 < \delta < 1$, we have for some constants $C_\pm(\alpha, u^0)$,

$$ \frac{C_-(\alpha, u^0)}{\delta} \le \int_0^\infty t^{-\alpha-\delta} U(t) dt, \quad \int_1^\infty t^{-\alpha-\delta} U(t) dt \le M \frac{C_+(\alpha)}{\delta}. \quad (18) $$

**Proof.** From (16), we obtain (17) since

$$ \int_0^t \frac{U(\tau)}{(1+t-\tau)^\alpha} d\tau = \int_0^\infty u^0(a) da + o(1), \quad t \gg 1. $$
The inequality (18) follows by integrating (16) with the weight $t^{-(1+\delta)}$ for $t > 1$. For the upper bound, due to (16) and the positivity of $u^0(a)$ in (7), we argue as follows:

$$ \int_{1}^{\infty} t^{-(1+\delta)} \int_{0}^{t} \frac{U(\tau)}{(1+t-\tau)^{\alpha}} d\tau dt < \int_{1}^{\infty} t^{-(1+\delta)} dt \int_{0}^{\infty} u^{0}(a) da \leq \frac{1}{\delta} \int_{0}^{\infty} u^{0}(a) da. $$

Using the Fubini theorem, we find with a further reduction,

$$ \int_{1}^{\infty} U(\tau) \int_{\tau}^{\infty} \frac{t^{-(1+\delta)}}{(1+t-\tau)^{\alpha}} dt d\tau \leq \frac{1}{\delta} \int_{0}^{\infty} u^{0}(a) da. $$

from which we conclude the upper bound by using that, when $\tau \geq 1$,

$$ \frac{t^{-(1+\delta)}}{(1+t-\tau)^{\alpha}} \geq t^{-1-\alpha-\delta}. $$

We may also write

$$ \int_{1}^{\infty} t^{-(1+\delta)} \int_{0}^{t} \frac{U(\tau)}{(1+t-\tau)^{\alpha}} d\tau dt = \int_{0}^{\infty} U(\tau) \int_{\tau}^{\infty} \frac{t^{-(1+\delta)}}{(1+t-\tau)^{\alpha}} d\tau d\tau. $$

From (16), the left hand side of the above equation can be bounded from below such that

$$ \int_{1}^{\infty} t^{-(1+\delta)} \left[ \int_{0}^{\infty} u^{0}(a) da - \int_{0}^{\infty} u^{0}(a) \frac{(1+a)^{\alpha}}{(2+a)^{\alpha}} da \right] dt \leq \int_{1}^{\infty} t^{-(1+\delta)} \int_{0}^{t} \frac{U(\tau)}{(1+t-\tau)^{\alpha}} d\tau dt, $$

from which we conclude the lower bound in (18).

**Remark 4** The interpretation is that the density, which is conserved, is spreading to infinity as for the heat equation in the whole space. The factor $d(a)$ gives less weight for large $a$ and thus does not see any mass.

# 4 sub-diffusion equation

## 4.1 Assumptions and main result

Next, we consider the scenario where particles undergo spatial jumps upon renewal while inside a trap located at position $y$, the duration time of particles follows an age-structured model. Upon exiting the trap, the transition probability of jumping to a new position $x$ is governed by a rate $\omega_{\epsilon}(x, y)$, which satisfies the following conditions: it is a measurable function such that

$$ \omega_{\epsilon}(x, y) \geq 0, \quad \int_{\mathbb{R}^d} \omega_{\epsilon}(x, y) dy = 1, \quad \text{and} \quad \omega_{\epsilon}(x, y) = \omega_{\epsilon}(y, x). \tag{19} $$

Consequently, the age-structured jump model can be formulated as:

$$ \left\{ \begin{array}{l} \varepsilon^{\beta} \partial_t u_{\varepsilon}(t, a, x) + \partial_a u_{\varepsilon}(t, a, x) + d(a)u_{\varepsilon}(t, a, x) = 0, \quad a \ge 0, \ x \in \mathbb{R}^d, \ t \ge 0, \\ U_{\varepsilon}(t, x) := u_{\varepsilon}(t, 0, x) = \int_{0}^{\infty} \int_{\mathbb{R}^d} d(a)u_{\varepsilon}(t, a, y)\omega_{\varepsilon}(x, y)dady, \\ u_{\varepsilon}(0, a, x) = u^0(a, x). \end{array} \right. \tag{20} $$
As before, we define

$$D(a) = \int_{0}^{a} d(a')da' \quad (e^{-D(a)} = \frac{1}{(1+a)^{\alpha}} \text{ in case of (9)}).$$

We assume that the initial data satisfies, for some constant $C^0$,

$$u^0 \ge 0, \quad e^{D(a)}u^0 \in L^1((0, \infty) \times \mathbb{R}^n), \quad u^0(a, x) \le C^0 e^{-D(a)}. \quad (21)$$

The parameter $\epsilon$ is defined from the scaled probability that a particles at position $y$ will jump to position $x$. We assume that, additionally to (19), the transition probability $\omega_{\epsilon}(x, y)$ satisfies, as $\epsilon \to 0$,

$$\left\{ \begin{array}{ll} \displaystyle \int_{\mathbb{R}^d} \frac{1}{\epsilon^2} (x-y) \otimes (x-y) \omega_{\epsilon}(x,y) dy \to A(x) > 0 & \text{a symmetric positive matrix,} \\ \\ \displaystyle \frac{1}{\epsilon^3} \int_{\mathbb{R}^d} |x-y|^3 \omega_{\epsilon}(x,y) dy \in L^{\infty}(\mathbb{R}^d) & \text{uniformly in } \epsilon. \end{array} \right. \quad (22)$$

In one dimension, a simple example of $\omega_{\epsilon}$ that satisfies the above requirements is

$$\omega_{\epsilon}(x, y) = \frac{1}{\epsilon} f_{\epsilon}(\frac{|x-y|}{\epsilon}) = \frac{1}{\epsilon} f_{\epsilon}(|z|) = \begin{cases} \frac{1-|z|}{\epsilon}, & z \in [-1, 1), \\ 0, & z \in (-\infty, -1) \cup (1, +\infty) \end{cases}$$

The value of $\beta$ in Eq. (20) depends on the jump rate $d(a)$ and is determined further.

With these assumptions, the problem is conservative and we fix the total mass $M$ defined as

$$\int_{0}^{\infty} \int_{\mathbb{R}^d} u_{\epsilon}(t, a, x) dadx = \int_{0}^{\infty} \int_{\mathbb{R}^d} u^{0}(a, x) dadx := M, \quad \forall t \ge 0.$$

Since $e^{-D(a)}$ is a steady state solution, the comparison principle shows that the solution satisfies,

$$u_{\epsilon}(t, a, x) \le C^0 e^{-D(a)}, \quad U_{\epsilon}(t, x) \le C^0, \quad \forall t \ge 0.$$

Also, with further assumptions on the initial data, further regularity is available, see subsection 4.3. Let us set

$$\rho_{\epsilon}(t, x) = \int_{0}^{\infty} u_{\epsilon}(t, a, x) da, \quad \rho^{0}(x) = \int_{0}^{\infty} u^{0}(a, x) da$$

and to simplify notations,

$$\nabla\nabla(A\rho) := \sum_{i,j} \partial_{i,j}^2 [A_{i,j}\rho].$$

We consider the case when $d(a)$ is determined by (9), i.e., $D(a) = \frac{1}{(1+a)^{\alpha}}$, or when $d(a) = D_0$, i.e., $D(a)e^{D_0 a}$ (this is simpler than taking $\alpha > 1$ and has the same effect). Consider the weak limit, in bounded measures, $\rho_0$ of $\rho_{\epsilon}$, as $\epsilon \to 0$ (a priori extraction of a subsequence is needed but uniqueness of the limit shows that the full family converges). The parameter $\alpha$ will decide of the time scale $\beta$ and of the limiting equation, namely, between normal diffusion when $d(a) = D_0$:

$$\left\{ \begin{array}{l} \partial_t \rho_0(t, x) - D_0 \nabla\nabla[A(x)\rho_0] = 0, \\ \rho_0(0, x) = \rho^0(x), \end{array} \right. \quad (23)$$
and sub-diffusion when $d(a)$ is determined by (9) with $0 < \alpha < 1$:

$$ \left\{ \begin{array}{l} \partial_t \int_0^t \frac{\rho_0(\tau, x)}{(t-\tau)^\alpha} d\tau - \nabla\nabla[A(x)\rho_0(t, x)] = t^{-\alpha}\rho^0(x), \\ \rho_0(0, x) = 0 \quad \text{for } t < 0. \end{array} \right. \qquad (24) $$

**Remark 5** It is interesting to note that the initial contribution of $u^0(a, x)$ appears solely in the source term rather than in the initial condition, which differs from the classical normal diffusion case. When $t \to 0+$, the dominant term on the left-hand side of (24) becomes $t^{-\alpha}\rho_0(0, x)$. Comparing with the right-hand side of (24), this indicates that $\rho_0(0+, x) = \rho^0(x)$, meaning $\rho_0(t, x)$ jumps from 0 to $\rho^0(x)$ at $t = 0$.

For the sub-diffusion case, we also have, when $\varepsilon \to 0$,

$$ \int_0^t \frac{U_0(\tau, x)}{(t-\tau)^\alpha} d\tau = U_0(0, x)O(t^{1-\alpha}) \quad \text{for } t \approx 0. $$

These equations are understood in the weak form, testing against a smooth test function $\psi(t, x)$ vanishing at $t = T$ and compactly supported in $\mathbb{R}^n$. For Eq. (23) this is

$$ \int_0^T \int_{\mathbb{R}^n} \rho_0(t, x) [-\partial_t \psi(t, x) - D_0 \sum_{i,j} A_{i,j}(x) \partial_{i,j}^2 \psi(t, x)] dxdt = \int_{\mathbb{R}^n} \psi(0, x) \rho_0(0, x) dx. $$

And for the sub-diffusion case, it is

$$ \int_0^T \int_{\mathbb{R}^n} \rho_0(t, x) [-\int_t^T \frac{\partial_\tau \psi(\tau, x)}{(\tau-t)^\alpha} d\tau - \sum_{i,j} A_{i,j}(x) \partial_{i,j}^2 \psi(t, x)] dxdt = \int_{t=0}^T \int_{a=0}^\infty \int_{\mathbb{R}^n} t^{-\alpha} \psi(t, x) u^0(a, x) dxdadt. $$

**Theorem 6 (Derivation of the macroscopic limit)** Assume (21), (19), (22) and choose $d(a)$ as in (9). Then, we have

(i) for $\alpha > 1$ and $\beta = 2$, the limit of $\rho_\varepsilon(t, x)$ satisfies the diffusion equation (23),

(ii) for $\alpha \in (0, 1)$ and $\alpha\beta = 2$, the limit of $\rho_\varepsilon(t, x)$ satisfies the sub-diffusion equation (24).

## 4.2 Proof of the main theorem

**Proof.** We argue with the Laplace transform of Eq. (20), which gives

$$ \partial_a \hat{u}_\varepsilon(s, a, x) + [\varepsilon^\beta s + d(a)] \hat{u}_\varepsilon(s, a, x) = \varepsilon^\beta u^0(a, x), $$

With $\hat{U}_\varepsilon(s, x) = \hat{u}_\varepsilon(s, 0, x)$ and the definition of $D(a)$, the solution is given by

$$ \hat{u}_{\varepsilon}(s, a, x) = \hat{U}_{\varepsilon}(s, x) e^{-D(a)-\varepsilon^{\beta}sa} + \varepsilon^{\beta} e^{-D(a)-\varepsilon^{\beta}sa} \int_{0}^{a} u^{0}(a', x) e^{D(a')+\varepsilon^{\beta}sa'} da'. \quad (25) $$

From this representation formula, we may deduce two macrocospic identities. Firstly, integrating in $a$, we find

$$ \hat{\rho}_{\varepsilon}(s, x) = \hat{U}_{\varepsilon}(s, x) \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da + \varepsilon^{\beta} \mathrm{ID}_{\varepsilon}^{(1)}(s, x) \quad (26) $$
where, changing $a$ to $a' + t$, we have

$$ 
\begin{aligned}
\mathrm{ID}_{\varepsilon}^{(1)}(s, x) &= \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} \int_{0}^{a} u^{0}(a', x) e^{D(a')+\varepsilon^{\beta}sa'} da' da \\
&= \int_{0}^{\infty} u^{0}(a', x) \int_{0}^{\infty} e^{D(a')-D(a'+t)-\varepsilon^{\beta}st} dt da'.
\end{aligned}
\quad (27) $$

Secondly, using the definition of $U_{\varepsilon}(t, x) = u_{\varepsilon}(t, 0, x)$ in (20), we multiply both sides of (25) (written in the variable $y$ rather than $x$) by $d(a)\omega_{\varepsilon}(x, y)$, integrate with respect to $a$ from $0$ to $\infty$ and $y$ in $\mathbb{R}^d$. This gives

$$ 
\begin{aligned}
\hat{U}_{\varepsilon}(s, x) &= \int_{0}^{\infty} \int_{\mathbb{R}^n} d(a) e^{-D(a)-\varepsilon^{\beta}sa} \hat{U}_{\varepsilon}(s, y) \omega_{\varepsilon}(x, y) dady + \varepsilon^{\beta} \mathrm{ID}_{\varepsilon}^{(2)}(s, x) \\
&= \int_{\mathbb{R}^n} \hat{U}_{\varepsilon}(s, y) \left[ 1 - \varepsilon^{\beta} s \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da \right] \omega_{\varepsilon}(x, y) dy + \varepsilon^{\beta} \mathrm{ID}_{\varepsilon}^{(2)}(s, x),
\end{aligned}
\quad (28) $$

where the initial data term is given by

$$ 
\begin{aligned}
\mathrm{ID}_{\varepsilon}^{(2)}(s, x) &= \int_{0}^{\infty} \int_{\mathbb{R}^n} d(a) e^{-D(a)-\varepsilon^{\beta}sa} \int_{0}^{a} u^{0}(a', y) e^{D(a')+\varepsilon^{\beta}sa'} da' \omega_{\varepsilon}(x, y) dady \\
&= \int_{a'=0}^{\infty} \int_{\mathbb{R}^n} u^{0}(a', y) e^{D(a')+\varepsilon^{\beta}sa'} \omega_{\varepsilon}(x, y) \int_{a'}^{\infty} d(a) e^{-D(a)-\varepsilon^{\beta}sa} dada' dy \\
&= \int_{a'=0}^{\infty} \int_{\mathbb{R}^n} u^{0}(a', y) \omega_{\varepsilon}(x, y) \left[ 1 - \varepsilon^{\beta} s \int_{a'}^{\infty} e^{-D(a)+D(a')-\varepsilon^{\beta}s(a-a')} \right] dada' dy.
\end{aligned}
\quad (29) $$

Our next step is to eliminate $\hat{U}_{\varepsilon}(s, x)$ from the above two formulas (26) and (28). Multiplying both sides of (26) by $\int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da$, we obtain,

$$ 
\begin{aligned}
\hat{\rho}_{\varepsilon}(s, x) - \varepsilon^{\beta} \mathrm{ID}_{\varepsilon}^{(1)}(s, x) &= \hat{U}_{\varepsilon}(s, x) \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da \\
&= \int_{\mathbb{R}^d} [\rho_{\varepsilon}(s, y) - \varepsilon^{\beta} \mathrm{ID}_{\varepsilon}^{(1)}(s, y)] [1 - \varepsilon^{\beta} s \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da] \omega_{\varepsilon}(x, y) dy \\
&\quad + \varepsilon^{\beta} \mathrm{ID}_{\varepsilon}^{(2)}(s, x) \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da.
\end{aligned}
$$

which we may rewrite as

$$ 
\begin{aligned}
\int_{\mathbb{R}^n} [\hat{\rho}_{\varepsilon}(s, x) - \hat{\rho}_{\varepsilon}(s, y) + \varepsilon^{\beta} s \hat{\rho}_{\varepsilon}(s, y) \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da] \omega_{\varepsilon}(x, y) dy \\
= Q_{\varepsilon}^{(2)}(s, x) + Q_{\varepsilon}^{(1)}(s, x).
\end{aligned}
\quad (30) $$

with

$$ 
\left\{
\begin{array}{l}
Q_{\varepsilon}^{(2)}(s, x) = \varepsilon^{\beta} \mathrm{ID}_{\varepsilon}^{(2)}(s, x) \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da, \\
Q_{\varepsilon}^{(1)}(s, x) = \varepsilon^{\beta} \int_{\mathbb{R}^d} \left[ \mathrm{ID}_{\varepsilon}^{(1)}(s, x) - \mathrm{ID}_{\varepsilon}^{(1)}(s, y) \left[ 1 - \varepsilon^{\beta} s \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da \right] \right] \omega_{\varepsilon}(x, y) dy.
\end{array}
\right.
\quad (31) $$

We now analyze the different terms arising in the identity (30), focusing on the case when $d(a)$ is given by Eq. (9). We begin with an observation.
**General calculations.** We first recall that

$$ 
\begin{aligned}
\varepsilon^{\beta} \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da &= \varepsilon^{\beta} \int_{0}^{\infty} \frac{e^{-\varepsilon^{\beta}sa}}{(1+a)^{\alpha}} da = \varepsilon^{\alpha\beta} \int_{0}^{\infty} \frac{e^{-st}}{(\varepsilon^{\beta}+t)^{\alpha}} dt \\
&= \varepsilon^{\alpha\beta} \int_{0}^{\infty} \frac{e^{-st}}{t^{\alpha}} dt [1+O(\varepsilon^{\beta})] = \varepsilon^{\alpha\beta} \mathcal{L}_s(t^{-\alpha}) [1+O(\varepsilon^{\beta})].
\end{aligned}
\quad (32) $$

In a similar way, with the change of variable $a \to t = \varepsilon^{\beta}(a - a')$

$$ 
\begin{aligned}
\varepsilon^{\beta} \int_{a'}^{\infty} e^{D(a')-D(a)-\varepsilon^{\beta}s(a-a')} da &= \varepsilon^{\alpha\beta} \int_{0}^{\infty} \frac{(1+a')^{\alpha}}{(t+(1+a')\varepsilon^{\beta})^{\alpha}} e^{-st} dt \\
&= \varepsilon^{\alpha\beta} \int_{0}^{\infty} \frac{(1+a')^{\alpha}}{t^{\alpha}} e^{-st} dt [1+O(\varepsilon^{\beta})] = \varepsilon^{\alpha\beta} (1+a')^{\alpha} \mathcal{L}_s(t^{-\alpha}) [1+O(\varepsilon^{\beta})].
\end{aligned}
\quad (33) $$

**The terms for $\hat{\rho}_{\varepsilon}$.** We begin with the first two terms on the left hand side (LHS) of (30). With assumptions (19), in the weak sense, by testing against a smooth test function with compact support $\psi(x)$, they give

$$ 
\begin{aligned}
\int_{\mathbb{R}^{2n}} \psi(x)[\hat{\rho}_{\varepsilon}(s, y) - \hat{\rho}_{\varepsilon}(s, x)]\omega_{\varepsilon}(x, y)dydx &= \int_{\mathbb{R}^{2n}} \hat{\rho}_{\varepsilon}(s, x)[\psi(y) - \psi(x)]\omega_{\varepsilon}(x, y)dydx \\
&= \varepsilon^2 \int_{\mathbb{R}^{2n}} \hat{\rho}_{\varepsilon}(s, x)[\partial_{i,j}^2\psi(x)(y_i - x_i)(y_j - x_j) + O(|x-y|^3)]\omega_{\varepsilon}(x, y)dydx \\
&= \varepsilon^2 \int_{\mathbb{R}^{n}} \hat{\rho}_{\varepsilon}(s, x)[A(x) : \nabla\psi(x) + O(\varepsilon)]dx.
\end{aligned}
$$

The third term on the (LHS) of (30), can be treated with the same argument and yields

$$ \varepsilon^{\beta} s \int_{\mathbb{R}^n} \hat{\rho}_{\varepsilon}(s, y) \int_0^\infty e^{-D(a)-\varepsilon^{\beta}sa} da \omega_{\varepsilon}(x, y) dy = \varepsilon^{\beta} s \hat{\rho}_{\varepsilon}(s, x) \int_0^\infty e^{-D(a)-\varepsilon^{\beta}sa} da [1+O(\varepsilon^2)] \quad (34) $$

Following the calculation for (12), when $d(a) = \frac{\alpha}{1+a}$, $0 < \alpha < 1$, the RHS of (34) can be written, with $\tau = \varepsilon^{\beta}a$ and according to (32)

$$ 
\begin{aligned}
\varepsilon^{\beta} s \int_{0}^{\infty} \rho_{\varepsilon}(x, t) e^{-st} dt \int_{0}^{\infty} \frac{e^{-\varepsilon^{\beta}sa}}{(1+a)^{\alpha}} da
&= s \int_{0}^{\infty} \rho_{\varepsilon} e^{-st} dt \int_{0}^{\infty} \frac{e^{-s\tau}}{(1+\varepsilon^{-\beta}\tau)^{\alpha}} d\tau \\
&= \varepsilon^{\alpha\beta} s \int_{0}^{\infty} \rho_{\varepsilon} e^{-st} \int_{t}^{\infty} \frac{e^{-s(\tau-t)}}{(\varepsilon^{\beta} + \tau - t)^{\alpha}} d\tau dt \\
&= \varepsilon^{\alpha\beta} s \int_{0}^{\infty} e^{-s\tau} \int_{0}^{\tau} \frac{\rho_{\varepsilon}(t)}{(\varepsilon^{\beta} + \tau - t)^{\alpha}} dt d\tau \\
&= \varepsilon^{\alpha\beta} \int_{0}^{\infty} e^{-st} \partial_t \int_{0}^{t} \frac{\rho_{\varepsilon}(\tau)}{(\varepsilon^{\beta} + t - \tau)^{\alpha}} d\tau dt.
\end{aligned}
$$

Here the last equality is due to integration by parts.

As a conclusion, we have obtained that, in the weak sense, the contribution of the left hand side of (30) is given by

$$ LHS = -\varepsilon^2 \nabla \nabla [A(x) \hat{\rho}_{\varepsilon}(s, x)] [1 + O(\varepsilon)] + \varepsilon^{\alpha\beta} \int_{0}^{\infty} e^{-st} \partial_t \int_{0}^{t} \frac{\rho_{\varepsilon}(\tau)}{(t-\tau)^{\alpha}} d\tau dt [1 + O(\varepsilon^{\beta})]. \quad (35) $$

At this stage, one can already observe that terms of orders $\varepsilon^2$ and $\varepsilon^{\alpha\beta}$ should be compatible and thus the choice $2 = \alpha\beta$.
**The term Q<sub>ε</sub><sup>(1)</sup>(s, x).** According to the analysis for ρ<sub>ε</sub>, the jump kernel ω<sub>ε</sub>(x, y) produces an O(ε<sup>2</sup>) compared to identity. Therefore, still in the weak sense, we have

$$
\begin{align*}
Q_{\varepsilon}^{(1)}(s,x) &= \varepsilon^{\beta} \left[ \mathrm{ID}_{\varepsilon}^{(1)}(s,x) - \mathrm{ID}_{\varepsilon}^{(1)}(s,x) [1 + O(\varepsilon^2)] [1 - \varepsilon^{\beta} s \int_{0}^{\infty} e^{-D(a)-\varepsilon^{\beta}sa} da] \right] \\
&= O(\varepsilon^{2+\beta}) \mathrm{ID}_{\varepsilon}^{(1)}(s,x) [1 + O(\varepsilon^{\alpha\beta})]
\end{align*}
$$

according to (32). The term ID<sub>ε</sub><sup>(1)</sup> can be also evaluated as

$$
\varepsilon^{\beta} \mathrm{ID}_{\varepsilon}^{(1)}(s,x) = \varepsilon^{\beta} \int_{0}^{\infty} \int_{0}^{\infty} u^{0}(a', x) \frac{(1+a')^{\alpha}}{(t+1+a')^{\alpha}} e^{-\varepsilon^{\beta}st} dt da' = O(\varepsilon^{\alpha\beta})
$$

according to (33).

As a consequence, the contribution of the term Q<sub>ε</sub><sup>(1)</sup>(s, x), vanishes compared to those arising in (35).

**The term Q<sub>ε</sub><sup>(2)</sup>(s, x).** To analyze the limit of Eq. (30), we finally need to analyze the term Q<sub>ε</sub><sup>(2)</sup>(s, x) in (31).

Using (32), we have

$$
Q_{\varepsilon}^{(2)}(s,x) = e^{\alpha\beta} \int_{0}^{\infty} \frac{e^{-st}}{t^{\alpha}} dt [1 + O(\varepsilon^{\beta})] \mathrm{ID}_{\varepsilon}^{(2)}(s,x).
$$

Also, according to the analysis for ρ<sub>ε</sub>, in the definition (29) of ID<sub>ε</sub><sup>(2)</sup>, the jump kernel ω<sub>ε</sub>(x, y) produces an O(ε<sup>2</sup>) compared to identity. Therefore, we also have

$$
\begin{align*}
\mathrm{ID}_{\varepsilon}^{(2)}(s,x) &= \int_{a'=0}^{\infty} u^0(a', x) [1 + O(\varepsilon^2)] [1 - \varepsilon^{\beta} s \int_{a'}^{\infty} e^{-D(a)+D(a')-\varepsilon^{\beta}s(a-a')} da'] da da' \\
&= \int_{0}^{\infty} u^0(a, x) (1+a)^{\alpha} [1 + O(\varepsilon^2)] [1 + O(\varepsilon^{\alpha\beta})] da,
\end{align*}
$$

where we have used the formula (33).

We conclude that

$$
Q_{\varepsilon}^{(2)}(s,x) = e^{\alpha\beta} \int_{0}^{\infty} \frac{e^{-st}}{t^{\alpha}} dt \int_{0}^{\infty} u^{0}(a,x) (1+a)^{\alpha} da [1 + O(\varepsilon^2 + \varepsilon^{\alpha\beta})]. \quad (36)
$$

**Conclusion (sub-diffusion case).** Back to Eq. (30), we choose αβ = 2. Putting together (35) and (36), we find

$$
\begin{align*}
& \int_0^\infty e^{-st} \partial_t \int_0^t \frac{\rho_\varepsilon(\tau, x)}{(t-\tau)^\alpha} d\tau dt \\
& = \nabla \nabla [A(x) \hat{\rho}_\varepsilon(s, x)] + \int_0^\infty t^{-\alpha} e^{-st} \int_0^\infty u^0(a, x) da dt + O(\varepsilon).
\end{align*}
$$

This is the Laplace transform of the equation

$$
\partial_t \int_0^t \frac{\rho_\epsilon(\tau, x)}{(t-\tau)^\alpha} d\tau = \nabla \nabla (A(x) \rho_\epsilon(s, x)) + t^{-\alpha} \int_0^\infty u^0(a, x) da + O(\epsilon). \quad (37)
$$

As ε vanishes, we find that the limit ρ₀ of ρ<sub>ε</sub> (which exists by uniqueness), satisfies Eq. (24).
**The diffusion case.** When $d(a) := D_0$, where $D_0$ is a positive constant, then $D(a) = D_0 a$. We can compute explicitly

$$ \int_0^\infty e^{-D_0 a - \varepsilon^\beta s a} da = \frac{1}{D_0 + \varepsilon^\beta s}, \quad \int_{a'}^\infty e^{-D(a)+D(a')-\varepsilon^\beta s(a-a')} da = \frac{1}{D_0 + \varepsilon^\beta s}, $$

when $d(a) = \frac{\alpha}{1+a}$ with $\alpha > 1$, similar formula hold with $e^{-D(a)} = \frac{1}{(1+a)^\alpha}$ and $D_0 = 1 - \alpha$ in the above equalities, at leading order. Then the expressions in (30) are simpler and reads at higher order

$$ \varepsilon^2 \nabla \nabla (A(x) \hat{\rho}_\varepsilon(s, x)) + \varepsilon^\beta \frac{s}{D_0} \hat{\rho}_\varepsilon(s, x) = Q_\varepsilon^{(1)}(s, x) + Q_\varepsilon^{(2)}(s, x), $$

from which we conclude that the scale is now $\beta = 2$.

The terms on the right hand side can be evaluated easily using their definitions in (31). For $Q_\varepsilon^{(1)}$, we find

$$ Q_\varepsilon^{(1)}(s, x) = \varepsilon^\beta O(\varepsilon^2) \mathrm{ID}_\varepsilon^{(1)}(s, x) $$

and, using the definition in (27), at leading order

$$ \mathrm{ID}_\varepsilon^{(1)}(s, x) = \frac{1}{D_0} \int_0^\infty u^0(a', x). $$

Therefore the contribution of $Q_\varepsilon^{(1)}(s, x)$ is again at smaller scale.

The term $Q_\varepsilon^{(2)}$ yields

$$ Q_\varepsilon^{(2)}(s, x) = \frac{\varepsilon^\beta}{D_0} \mathrm{ID}_\varepsilon^{(2)}(s, x) = \frac{\varepsilon^\beta}{D_0} \int_0^\infty u^0(a, x) da [1 + O(\varepsilon^\beta)]. $$

Finally, we find, in the weak sense

$$ \varepsilon^2 \nabla \nabla (A(x) \hat{\rho}_\varepsilon(s, x)) + \varepsilon^\beta \frac{s}{D_0} \hat{\rho}_\varepsilon(s, x) = \frac{\varepsilon^\beta}{D_0} \int_0^\infty u^0(a, x) da + O(\varepsilon^{2\beta}), $$

which is the Laplace transform of the equation

$$ \partial_t \rho_\varepsilon(t, x) = D_0 \nabla \nabla (A(x) \rho_\varepsilon(t, x)) + \delta(t) \rho^0(x) + O(\varepsilon^2) $$

Here we recall our convention that $u(t, a, x)$ and $\rho(t, x)$ are extended by 0 for $t < 0$ so that the Dirac mass represents the initial data as a jump. Consequently, the quantity $\rho^0(x) = \int_0^\infty u^0(x, a) da$ represents the initial condition for $\rho_0(t, x)$.

## 4.3 Remarks on regularity

**Regularity of $\rho_0$ in the sub-diffusion case.** When $\alpha < 1$, the only uniform bound available for $\rho_\varepsilon$ is that $\rho_\varepsilon \ge 0$ and $\int_{\mathbb{R}^d} \rho_\varepsilon(t, x) dx = M$. Consequently the limit $\rho_0$ is a measure. However it is standard to obtain some Sobolev estimates in $L^2$ (see [9] for further regularity). We consider the case when $A = \mathrm{Id}$, the identity matrix, to simplify and explain now how to find the following $H^1$ bound for $\rho_0$:

$$ \begin{aligned} & \frac{1}{2} \partial_t \int_0^t \int_{\mathbb{R}^d} \frac{v^2(\tau, x)}{(t-\tau)^\alpha} d\tau dx + \int_{\mathbb{R}^d} |\nabla v|^2 dx + \frac{\alpha}{2} \int_0^t \int_{\mathbb{R}^d} \frac{|v(\tau, x) - v(t, x)|^2}{(t-\tau)^{\alpha+1}} d\tau dx + \int_{\mathbb{R}^d} \frac{v^2(t)}{2t^\alpha} dx \\ & = - \int_{\mathbb{R}^d} \nabla v(t, x) \cdot \nabla \rho^0(x) dx, \end{aligned} \quad (38) $$
where $v(t, x) = \rho_0(t, x) - \rho^0(x)$.

Notice that $v$ satisfies the sub-diffusion equation

$$ \partial_t \int_0^t \frac{v(\tau)}{(t-\tau)^\alpha} d\tau - \Delta v(t, x) = \Delta \rho^0(x). $$

From Remark 5, $v(0, x) = 0$. Then, (38) is obtained by multiplying this equation by $v(t, x)$ and using integration by parts. The fractional time derivative term follows from the following identity: for $v \in C^1([0, \infty))$ satisfying $v(0) = 0$, one has

$$ \frac{1}{2} \partial_t \int_0^t \frac{v^2(\tau)}{(t-\tau)^\alpha} d\tau + \frac{v^2(t)}{2t^\alpha} + \frac{\alpha}{2} \int_0^t \frac{|v(\tau) - v(t)|^2}{(t-\tau)^{\alpha+1}} d\tau = v(t) \partial_t \int_0^t \frac{v(\tau)}{(t-\tau)^\alpha} d\tau. \quad (39) $$

To show this, on the one hand, we write

$$ \begin{aligned} \partial_t \int_0^t \frac{v(\tau)}{(t-\tau)^\alpha} d\tau &= \partial_t \int_0^t \frac{v(t-s)}{s^\alpha} ds = - \int_0^t \frac{\partial_s [v(t-s) - v(t)]}{s^\alpha} ds \\ &= - \left. \frac{v(t-s) - v(t)}{s^\alpha} \right|_{s=0}^{s=t} - \alpha \int_0^t \frac{v(t-s) - v(t)}{s^{\alpha+1}} ds \\ &= \frac{v(t)}{t^\alpha} + \lim_{s \to 0} \frac{v(t-s) - v(t)}{s} s^{1-\alpha} - \alpha \int_0^t \frac{v(t-s) - v(t)}{s^{\alpha+1}} ds \\ &= \frac{v(t)}{t^\alpha} - \alpha \int_0^t \frac{v(t-s) - v(t)}{s^{\alpha+1}} ds. \end{aligned} \quad (40) $$

On the other hand, applying the result to $v^2(\tau)$ we write

$$ \begin{aligned} \partial_t \int_0^t \frac{v(\tau)^2}{(t-\tau)^\alpha} d\tau &= \frac{v(t)^2}{t^\alpha} - \alpha \int_0^t \frac{v(t-s)^2 - v(t)^2}{s^{\alpha+1}} ds \\ &= \frac{v(t)^2}{t^\alpha} - \alpha \int_0^t \frac{|v(t-s) - v(t)|^2}{s^{\alpha+1}} ds - 2v(t)\alpha \int_0^t \frac{v(t-s) - v(t)}{s^{\alpha+1}} ds. \end{aligned} $$

Thanks to (40), we conclude

$$ \partial_t \int_0^t \frac{v^2(\tau)}{(t-\tau)^\alpha} d\tau = \frac{v^2(t)}{t^\alpha} - \alpha \int_0^t \frac{|v(t-s) - v(t)|^2}{s^{\alpha+1}} ds + 2v(t) \left[ \partial_t \int_0^t \frac{v(\tau)}{(t-\tau)^\alpha} d\tau - \frac{v(t)}{t^\alpha} \right], $$

which gives the announced result Eq. (39).

As in [1], notice that it also holds, for all convex function $\Phi$ such that $\Phi(0) = 0$,

$$ \partial_t \int_0^t \frac{\Phi(v(\tau))}{(t-\tau)^\alpha} d\tau \le \Phi'(v(t)) \partial_t \int_0^t \frac{v(\tau)}{(t-\tau)^\alpha} d\tau. $$

It would be interesting to understand if a related estimate holds uniformly in $\varepsilon$.

**Time regularity.** The time regularity for $U(t, x)$ can also be obtained with the following observation. Let $z_\varepsilon = \partial_t u_\varepsilon$. Differentiating in time on both sides of Eq. (20), we find that $z_\varepsilon$ satisfies the same equation. Therefore, we also obtain

$$ |\partial_t z_\varepsilon(t, a, x)| \le \frac{C^1}{(1+a)^\alpha} $$
under the assumption $|\partial_t z_\epsilon(0, a, x)| \le \frac{C^1}{(1+a)^\alpha}$, which, using Eq. (20), requires that

$$ \partial_a u_\epsilon^0 + d(a)u_\epsilon^0 = O(\epsilon^\beta). \qquad (41) $$

This 'near equilibrium' condition implies, after integration in y and a against the weight $\frac{\omega_\epsilon(x,y)}{1+a}$, that the renewal term $U_\epsilon(t, x)$ satisfies

$$ |\partial_t U_\epsilon(t, x)| \le \frac{C^1}{\alpha}. $$

## References

[1] M. ALLEN, L. CAFFARELLI, AND A. VASSEUR, Porous medium flow with both a fractional potential pressure and fractional time derivative, Chin. Ann. Math., 38 (2017), p. 45–82.

[2] C. N. ANGSTMANN, I. C. DONNELLY, B. I. HENRY, T. A. M. LANGLANDS, AND P. STRAKA, Generalized continuous time random walks, master equations, and fractional Fokker-Planck equations, SIAM Journal on Applied Mathematics, 75 (2015), pp. 1445–1468.

[3] E. BARKAI, R. METZLER, AND J. KLAFTER, From continuous time random walks to the fractional Fokker-Planck equation, Phys. Rev. E, 61 (2000), pp. 132–138.

[4] H. BERRY, T. LEPOUTRE, AND A. M. GONZÁLEZ, Quantitative convergence towards a self-similar profile in an age-structured renewal equation for subdiffusion, Acta Appl. Math., 145 (2016), pp. 15–45.

[5] E. CAGLIOTI AND F. GOLSE, On the Boltzmann-Grad limit for the two-dimensional periodic Lorentz gas, J. Stat. Phys., 141 (2010), pp. 264–317.

[6] V. CALVEZ, P. GABRIEL, AND Á. M. GONZÁLEZ, Limiting Hamilton-Jacobi equation for the large scale asymptotics of a subdiffusion jump-renewal equation, Asymptotic Anal., 115 (2019), pp. 63–94.

[7] M. CAPUTO, Diffusion of fluids in porous media with memory, Geothermics, 28 (1999), p. 113–130.

[8] O. DIEKMANN AND J. A. P. HEESTERBEEK, Mathematical epidemiology of infectious diseases, Wiley Series in Mathematical and Computational Biology, John Wiley & Sons, Ltd., Chichester, 2000. Model building, analysis and interpretation.

[9] H. DONG AND D. KIM, $L_p$-estimates for time fractional parabolic equations in divergence form with measurable coefficients, J. Funct. Anal., 278 (2020), p. 66. Id/No 108338.

[10] C. DÜLL, P. GWIAZDA, A. MARCINIAK-CZOCHRA, AND J. SKRZECZKOWSKI, Spaces of measures and their applications to structured population models, vol. 36 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, 2022.

[11] G. ESTRADA-RODRIGUEZ, H. GIMPERLEIN, K. J. PAINTER, AND J. STOCEK, Space-time fractional diffusion in cell movement models with delay, Math. Models Methods Appl. Sci., 29 (2019), pp. 65–88.
[12] M. FRANK AND T. GOUDON, *On a generalized Boltzmann equation for non-classical particle transport*, Kinet. Relat. Models, 3 (2010), pp. 395–407.

[13] M. FRANK AND W. SUN, *Fractional diffusion limits of non-classical transport equations*, Kinet. Relat. Models, 11 (2018), pp. 1503–1526.

[14] M. FRITZ, E. SÜLI, AND B. WOHLMUTH, *Analysis of a dilute polymer model with a time-fractional derivative*, SIAM J. Math. Anal., 56 (2024), pp. 2063–2089.

[15] M. IANNELLI AND F. MILNER, *The basic approach to age-structured population dynamics*, Lecture Notes on Mathematical Modelling in the Life Sciences, Springer, Dordrecht, 2017. Models, methods and numerics.

[16] J.-H. JEON, V. TEJEDOR, S. BUROV, E. BARKAI, C. SELHUBER-UNKEL, K. BERG-SØRENSEN, L. ODDERSHEDE, AND R. METZLER, *In vivo anomalous diffusion and weak ergodicity breaking of lipid granules*, Physical review letters, 106 (2011), p. 048103.

[17] J. MARKLOF AND A. STRÖMBERGSSON, *The distribution of free path lengths in the periodic Lorentz gas and related lattice point problems*, Ann. Math. (2), 172 (2010), pp. 1949–2033.

[18] Y. MEROZ AND I. M. SOKOLOV, *A toolbox for determining subdiffusive mechanisms*, Physics Reports, 573 (2015), pp. 1–29.

[19] R. METZLER AND J. KLAFTER, *The random walk's guide to anomalous diffusion: a fractional dynamics approach*, Physics Reports, 339 (2000), pp. 1–77.

[20] B. PERTHAME, *Transport equations in biology*, Frontiers in Mathematics, Birkhäuser Verlag, Basel, 2007.

[21] S. A. TABEI, S. BUROV, H. Y. KIM, A. KUZNETSOV, T. HUYNH, J. JURELLER, L. H. PHILIPSON, A. R. DINNER, AND N. F. SCHERER, *Intracellular transport of insulin granules is a subordinated random walk*, Proceedings of the National Academy of Sciences, 110 (2013), pp. 4911–4916.

[22] P. TAN, Y. LIANG, Q. XU, E. MAMONTOV, J. LI, X. XING, AND L. HONG, *Gradual crossover from subdiffusion to normal diffusion: A many-body effect in protein surface water*, Phys. Rev. Lett., (2018), p. 248101.

[23] ——, *Many-body effect renders universal subdiffusion to water on different proteins*, Phys. Rev. Research, (2022), p. L022003.

[24] H. R. THIEME, *Mathematics in population biology*, Princeton Series in Theoretical and Computational Biology, Princeton University Press, Princeton, NJ, 2003.

[25] I. M. TOLIĆ-NØRRELYKKE, E.-L. MUNTEANU, G. THON, L. ODDERSHEDE, AND K. BERG-SØRENSEN, *Anomalous diffusion in living yeast cells*, Physical review letters, 93 (2004), p. 078102.

[26] M. WEISS, M. ELSNER, F. KARTBERG, AND T. NILSSON, *Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells*, Biophysical journal, 87 (2004), pp. 3518–3524.

[27] R. ZACHER, *Global strong solvability of a quasilinear subdiffusion problem*, J. Evol. Equ., 12 (2012), pp. 813–831.