Subdiffusion in the Presence of Reactive Boundaries: A Generalized Feynman–Kac Approach

We derive, through subordination techniques, a generalized Feynman–Kac equation in the form of a time fractional Schrödinger equation. We relate such equation to a functional which we name the subordinated local time. We demonstrate through a stochastic treatment how this generalized Feynman–Kac equation describes subdiffusive processes with reactions. In this interpretation, the subordinated local time represents the number of times a specific spatial point is reached, with the amount of time spent there being immaterial. This distinction provides a practical advance due to the potential long waiting time nature of subdiffusive processes. The subordinated local time is used to formulate a probabilistic understanding of subdiffusion with reactions, leading to the well known radiation boundary condition. We demonstrate the equivalence between the generalized Feynman–Kac equation with a reflecting boundary and the fractional diffusion equation with a radiation boundary. We solve the former and find the first-reaction probability density in analytic form in the time domain, in terms of the Wright function. We are also able to find the survival probability and subordinated local time density analytically. These results are validated by stochastic simulations that use the subordinated local time description of subdiffusion in the presence of reactions.


Introduction
In recent years anomalous diffusion has been found to be a prevalent transport mechanism across many systems.Specifically, subdiffusion is of key interest due to its defining sub-linear mean-square displacement in time, i.e.
where Y (t) is a time dependent random variable with α ∈ (0, 1).Due to this sub-linear form, subdiffusive motion has been observed in a variety of physical and biological processes (see [1,2] and references therein).Over the last two decades or so, much work has been endeavored to create a unified framework to describe subdiffusive motion.One of the most utilised approach is a fractional diffusion or Fokker-Planck equation, derived from a generalized master equation (GME) or continuous-time random walk (CTRW) approach [3][4][5].It is also possible to obtain this fractional diffusion equation through a subordinated Langevin approach [6].
Within this unified framework it is natural to consider another fundamental equation in the study of stochastic processes, the Feynman-Kac equation (FKE), and its extension to the subdiffusive case.The classical FKE is a well known tool to study functionals of Brownian motion with numerous applications across physics and other areas of science [7].Thus there was a clear need for the extension to when the underlying stochastic process is subdiffusive.This need has been met in recent years with a fractional FKE having been found to study functionals of subdiffusion [8][9][10][11] with further generalizations to space and time dependent forces [12,13], tempered subdiffusion [14], aging subdiffusion [15], multiplicative noise [16] and reaction-subdiffusion processes [17].
One of the most utilized functionals is the so-called local time functional [18], which finds applications in various areas [19], and has recently been used to build a probabilistic description of diffusion with surface reactions [20][21][22].In this approach surface reactions are described via stopping conditions based on the local time of the Brownian particle at the boundary.The Brownian particle undergoes normal diffusion being reflected each time it reaches the boundary until the local time at the boundary exceeds a random variable drawn from an exponential distribution with the inverse scale being the reactivity parameter.The time at which this occurs is then the reaction time, where the particle has then reacted, absorbed, changed species etc.This approach presents a formal and practical advance compared to classical methods such as using radiation boundary conditions [23] or placing partial traps or defects in the domain [24,25].
Due to the ubiquity of subdiffusion in complex systems, these systems are often bounded by reactive boundaries [26].Thus there is a clear need to generalize this description for when the motion is subdiffusive.The main purpose of this paper is to provide such a generalization.We do this by considering an alternative generalized FKE [27,28] rather than the previously mentioned fractional FKE.The generalized FKE we use is in the form of an (imaginary time) time fractional Schrödinger equation [29][30][31][32] and governs subordinated forms of the functionals.This proves to be a useful recipe in the case of the local time functional for providing such a generalized description of subdiffusion in the presence of reactive boundaries.
The paper is structured as follows.In Sec. 2 we recall the classical FKE to which we derive a generalized form through subordination techniques and introduce the subordinated local time functional whose meaning is uncovered using a CTRW approach.In Sec. 3 we present a probabilistic interpretation of this generalized FKE as subdiffusion with reactions using the subordinated local time and how this is connected to the radiation boundary condition (BC).In Sec. 4 we present an application of these findings by analytically studying three important quantities associated with subdiffusion in the presence of a radiation boundary, namely the first-reaction time density, survival probability and subordinated local time density.We confirm these analytic results with stochastic simulations.Finally, we discuss and conclude our findings in Sec. 5.

The Classical Feynman-Kac Equation
The celebrated FKE, derived in 1949 by Kac influenced by Feynman's path integral description of quantum mechanics, has become a fundamental tool in the theory of stochastic processes [33,34].The Feynman-Kac theory provides a rigorous connection between the paths, X(t), of a Bronwnian motion process and the solution to the (imaginary time) Schrödinger equation [35].The main utility however is the connection to functionals of Brownian motion [7], where U (x) is some arbitrary function.The FKE governs the (Laplace/Fourier transformed) joint probability density, ρ(x, A, t|x 0 ), of X(t) and A(t), given by [35,36] where K is the diffusion coefficient, i.e. it is the strength of the delta correlated noise of the Langevin equation associated with X(t), while the Laplace variable p is related to A via [7] It should be noted that if A(t) is not always positive, then the Laplace transform needs to be replaced by a Fourier transform, i.e. p → −ip and the lower integration bound changed to −∞ [9].Alternatively, Eq. ( 3) can be represented via the expectation, In Eq. ( 4) the average is over all trajectory realizations of X(t) that starts at X(0) = x 0 , that is P (x, p, 0|x 0 ) = δ(x − x 0 ).

Time-Changed Process
Let us now consider subdiffusion through a CTRW paradigm [37,38].A CTRW formalism is constructed by considering a random walker which waits at each step, i, for a time η i and then proceeds to jump a distance ξ i .The random variables ξ i and η i , are independent and identically distributed.Thus, after n steps, the position of the random walker, Y n and the total time elapsed, T n , can be found by [13], where Y 0 is the initial position.Through a parameterization of the CTRW via the continuous time variable, t, instead of the number of steps, n, Eq. ( 5) can be written compactly as, where N (t) = max {n ≥ 0 : T n ≤ t}, such that N (t) is a random variable itself, as a consequence of containing the statistics of the random waiting times.If we now take the continuum limit of Eq. ( 5), we obtain [13] X Here, τ is not the real physical time, but is instead an operational time.
In the same sense, we take the continuum limit of Eq. ( 6), with N (t) → S(t), to obtain Since S(0) = 0, we have In other words Y (t) has undergone a time change and is a subordinated process, such that S(t) can be interpreted as a stochastic clock [6].Specifically, subdiffusion is generated in the macroscopic limit of a CTRW with a distribution of waiting times that is heavy-tailed, such that the mean waiting time is infinite.If we indicate with T α (τ ) a waiting time distribution which follows an α-stble Lévy distribution, the Laplace transform of T α (τ ) is exp{−κT α (τ )} = exp{−τ κ α }, with α ∈ (0, 1) [39].The stochastic clock is then defined as, S α (t) = inf{τ > 0 : and will be termed the inverse α-stable subordinator [40][41][42].The Laplace transform of the probability density of S α (t), g(τ, t), is given by [43] g alternatively by taking the derivative of both sides of Eq. ( 10) and performing the inverse Laplace transform, one can show g(τ, t) satisfies the following fractional differential equation [43], Here 0 D 1−α t is a fractional derivative of Riemann-Liouville type [44], i.e. for a generic function f (t), When ξ(τ ) is the standard Langevin force (i.e.Gaussian white noise), Y (t) = X(S α (t)) is a subdiffusion process, with X(τ ) being standard Brownian motion with probability density, W (x, τ |x 0 ) = P (x, 0, τ |x 0 ).In Fig. (1) we show a simple realization of T α (τ ) and its corresponding S α (t), and how the resulting X(τ ) trajectory is modified to a Y (t) trajectory.
The probability density of Y (t) written as W α (x, t|x 0 ) can be given in terms of W (x, τ |x 0 ) and g(τ, t) as follows [46], due to the independence of X(τ ) and S α (t).By taking the time derivative of both sides of Eq. ( 13), whilst using Eq. ( 11) and integrating by parts we find, (10) it is clear that g(∞, t) = 0 and 0 D 1−α t g(0, t) = δ(t), thus for t > 0, and using the normal diffusion equation, we recover the fractional diffusion equation (FDE) [1,2,47], Fig. 1 A set of trajectories for the four process that leads to subdiffusion, (a) α-stable Lévy motion, (b) Brownian motion, (c) inverse α-stable subordinator, (d) subdiffusion.(a) and (b) are both in terms of the operational time τ and (c) and (d) are in terms of the real physical time t.These plots were generated using the algorithm in [45], where Tα(τ ) is simulated to which Sα(t) is found using Eq.(9).Then Y (t) = X(Sα(t)) can be approximated by interpolating X(τ ) for the values of Sα(t).One can see that the delays in (c) correspond to periods of waiting in (d).
with K α being the generalized diffusion coefficient which has dimensions [length] 2 /[time] α .Now let us consider not only subordinating the Brownian motion but also the functional of Brownian motion, i.e.

A(S
The joint density, P α (x, p, t|x 0 ), where will then be governed by a generalized FKE.Using the same arguments as above we have [27], It is then simple to find P α (x, p, t|x 0 ) via Eq.( 16), since Using Eq. ( 4) and the properties of independence, we have where ρ α (x, A, t|x 0 ) is the joint density of X(S α (t)) and A(S α (t)).
We point out the difference of Eq. ( 17) compared to the fractional FKE, as mentioned in Sec. 1, which is given by [8,9], where D 1−α t is the so-called fractional substantial derivative [48], Thus the corresponding density, ρ α (x, A, t|x 0 ), is then the joint density of Y (t) and A α (t) = t 0 U [Y (t )]dt (compare to Eq. ( 15)).Thus the functional is of the subdiffusive motion and is not a subordinated quantity, illustrating that Eq. ( 20) is the natural generalization of the FKE for studying functionals of subdiffusion.

Subordinated Local Time
The local time of a stochastic process, originally introduced by Lévy [18], is an important quantity that characterises the fraction of time the process spends at a certain point, x b .We will label the local time, (t), which is defined as the functional Without loss of generality, moving forward, we take this point to be at the origin, x b = 0. Now, let us consider the subordinated local time, (S α (t)), we are able to evince the meaning of this quantity, as follows.Returning to the CTRW paradigm (for the subdiffusive case), where the random walker is described by Eq. ( 6) we introduce the quantity N (t), which is the number of times the walker visits a region around the origin, ∂Ω, of width ε (ε gives the scale of the jump length i.e. ξ 2 i = ε 2 ), up to time t [49]: In Eq. ( 23) Y i is given by Eq. ( 5) and I ∂Ω (Y i ) is the indicator function, i.e I ∂Ω (Y i ) = 1 if Y i ∈ ∂Ω and 0 otherwise.Now we take a continuum limit such that Y n → X(τ ) and N (t) → S α (t) and introduce a scaling of K α /ε 2 to make N (t) dimensionless, Let us then take the diffusive limit, which entails letting ε vanish, resulting in the indicator function, I ∂Ω (X(s)), becoming the Dirac-δ function, then we have lim The limit in Eq. ( 25) exists due to the recurrent nature of Brownian motion in one dimension such that the number of visits, N (t), diverges in the diffusive limit and can be understood as the continuous analog of the (scaled) number of times a subdiffusive particle visits a certain point.Note this limit will not hold in higher dimensions, so one would need to consider a bounded domain where ∂Ω becomes a thin layer at a reflecting boundary [20,22].It is well known that the local time for a normal diffusive particle is the continuous limit of the number of times the particle visits a certain point [49].However, due to the long waiting times embedded in the subdiffusive dynamics a particle may spend anomalously long times at a certain point.Thus, the local time for a subdiffusive particle, α (t) = t 0 δ(Y (t ))dt , does not correspond to the continuous limit of the number of times the particle visits a point.So, if one only cares about whether the particle has reached a certain point a number of times, and not how long it has spent there, e.g. if an event occurs when the particle reaches that point for every visit, then the quantity of interest is the subordinated local time, Eq. ( 25), i.e. (S α (t)) = , so contrary to its name (S α (t)) is not actually a time.

Generalized Feynman-Kac Equation and Reactions
It is well known that the FKE (2) can be interpreted as diffusion with killing or reactions [36], such that the diffusive particle is removed from the system i.e. by being absorbed.Let us say for simplicity that absorption occurs at the origin, thus we take U (x) = δ(x).We introduce the killing or first reaction time T and assume the reaction dynamics are subordinated and thus governed by the stochastic clock, S α (t).Then the probability for a subdiffusive particle starting at x 0 to be killed in the time interval [t, t + h) can be approximated by Eq. ( 26) can be understood as follows.Every time the particle reaches the origin there is a chance that a reaction occurs and the trajectory will be killed.
As we are considering the reaction dynamics are occurring according to S α (t), the physical time interval [t, t+h) corresponds to the interval [S α (t), S α (t+h)) for the reaction dynamics.So for h << 1 one can approximate P[T ∈ [t, t+h)] x0 by multiplying the probability the particle is found at the origin at time t having not reacted previously, δ(X(S α (t))) x0 , by the interval over which the reactions occur, (S α (t + h) − S α (t)), and by the Laplace variable in Eq. ( 19), p, which here is considered as the reactivity.
As each reaction event is taken to be independent of each other, the probability the particle has not reacted (survived) up to time t is given by, where we have partitioned the interval [0, t] into 0 = t 0 < t 1 < ... < t M = t, where h = t i+1 − t i .In the limit we obtain, Due to the independence between the subdiffusive dynamics and the reaction events, we find that Eq. ( 19) (for U (x) = δ(x)), can be understood as [22, 50] This shows P α (x, p, t|x 0 ) can be interpreted as the probability density of a subdiffusive particle starting at x 0 to be at a position x whilst having not reacted at the origin.To link Eq. ( 29) to the joint probability density of the subordinated local time and position, ρ α (x, l, t|x 0 ), we use the integral form of P α (x, p, t|x 0 ) in Eq. ( 19) and with we integrate by parts to find Inside the integral we have the probability density of an exponentially distributed random variable, l, with mean 1/p, i.e.P[ l ∈ [l, l + dl)] = pe −pl dl.So if we replace pe −pl with δ(l − l) , we obtain Therefore, we have shown that the solution of the generalized FKE (17), P α (x, p, t|x 0 ), is the probability density of a subdiffusing particle to be found at a position x, whilst the subordinated local time has not exceeded the value of an exponentially distributed random variable l (see discussion in Sec. 1 and Refs.[20][21][22]50]).As (S α (t)) is a monotonically non-decreasing process, the event of { l > (S α (t))} is equivalent to {T > t}.Thus, by comparing Eq. ( 31) to Eq. ( 29) we can see how the reaction time, T , is related to the subordinated local time, (S α (t)), via In other words the reaction time and subordinated local time are intimately linked processes, with the reaction time being determined by the subordinated local time exceeding a certain value.

Radiation Boundary Condition
In the previous section we have established the connection between reactions and the generalized FKE, here we extend the idea to when the reactions are occurring on a boundary.We consider a CTRW generated by a nearestneighbour random walk moving on a discrete one-dimensional lattice with sites 0, 1, 2, ... and lattice spacing ε with the waiting time distribution, ψ(t), being heavy-tailed, thus ψ(t) υ α /t 1+α , where υ is a temporal scale parameter.The dynamics of the walker may be described by the GME for the occupation probability at lattice site i, initially starting from site j, W i,j (t), for i > 0 the GME is [51] where the Laplace transform of the memory kernel being, Φ( ) = ψ( )/(1 − ψ( )), here is again the Laplace variable.At the site i = 0, we have In the presence of a reactive boundary an incident particle at the lattice site i = 0, has a probability of reacting 1 − k and a probability of being reflected k.This situation can be summarized by the following flux condition, From Eq. ( 34) we can identify the discrete fluxes in and out of the boundary as J + 0,j (t) = t 0 dt Φ(t − t )W 0,j (t )/2 and J − 0,j (t) = t 0 dt Φ(t − t )W 1,j (t )/2, respectively [52].If we insert these fluxes into Eq.( 35), and using the following relation for the total flux, J 0,j (t) = J + 0,j (t) − J − 0,j (t), we obtain and find the relation in terms of W 0,j (t), as Let us now take the diffusive limit, which entails taking the limits ε → 0 and υ α → 0 [4].This corresponds to ψ( ) ∼ 1 − (υ ) α , so Φ( ) ∼ υ −α 1−α with W i,j (t)/ε → W α (x, t|x 0 ) and J i,j (t) → J α (x, t|x 0 ), where iε → x and jε → x 0 .By taking the Laplace transform of Eq. ( 37) and inserting Φ( ), we obtain where the flux is defined as x, t|x 0 ) (from writing Eq. ( 14) as a continuity equation, with ε 2 /(2υ α ) → K α .Eq. ( 39) implies that in the diffusive limit one requires k → 1 [53], which is a consequence of the number of visits to the origin becoming infinite.
Taking the inverse Laplace transform of Eq. ( 38), we have Clearly Eq. ( 40) is equivalent to the so-called radiation BC [23,26,38], The radiation BC (41) describes a reactive boundary such that an incident particle is either absorbed or reflected depending on the reactivity parameter λ α which has dimensions [length]/[time] α , such that λ α = ∞ represents full absorption and λ α = 0 represents full reflection.
Note that here we have obtained the reactivity parameter, λ α , derived from a reaction probability principle, whereas previously it had been found using a reaction rate [52,[54][55][56].

Connection Between the Generalized Feynman-Kac Equation and the Radiation Boundary Condition
Starting from the limit form of λ α in Eq. ( 39), for small ε the probability of reflection can be expressed as k ≈ 1/(1 + ελ α /K α ).The probability of visiting the origin n times without a reaction (being reflected) is simply k n since each interaction is independent, but with reference to Sec. 2.3, we know the particle visits the origin a random number of times, N (t).To find the probability of the particle, which started at x 0 , to have not reacted, we must average over all possible realisations of N (t), to which we obtain [57] P After inserting the expression for k into Eq.( 42) and expanding for small ε, we have Using Eq. ( 25), Eq. ( 43) in the limit ε → 0 becomes [57], which is exactly Eq. ( 28) with λ α in place of p.Note that because of the boundary X(τ ) is now reflected Brownian motion, whereas in Sec.3.1 we did not impose such a restriction.This implies that the solution to the FDE, (14) with the radiation BC ( 41) is equivalent to solving the generalized FKE (17) with U (x) = δ(x) and the reflecting (Neumann) BC, lim x→0 ∂ x P α (x, p, t|x 0 ) = 0. We can simply verify this by considering the generalized FKE with U (x) = δ(x), where we have replaced p with λ α as we are looking specifically at the radiation BC.Eq. ( 45) has been considered before in the specific context of geminate recombination [54,55,58].Integrating both sides of Eq. ( 45) over the range [0, ∆] with respect to x and making use of the reflecting BC, gives [55] ∆ 0 −λ α P α (0, λ α , t|x 0 ) .
By taking ∆ → 0 the left-hand side vanishes and the radiation BC ( 41) is recovered, showing the equivalent formulations of the problem.We note that it has recently been shown that another equation (not in the FKE form), coming from (normal) diffusion through permeable barriers, satisfies the radiation boundary condition [59] with certain conditions.We extend this equation to the subdiffusive case in Appendix A and show how it leads to the radiation BC.This fact further compounds the notion that the reaction dynamics is a subordinated process.
4 First-Reaction Time

Generalized Feynman-Kac Equation Solution
To study the first-reaction time (FRT) of a subdiffusive particle in the presence of a radiation boundary at the origin, we must find the corresponding probability density, P α (x, λ α , t|x 0 ).As demonstrated in Sec.3.3 this may be achieved through two equivalent methods, solving the FDE ( 14) with the radiation BC (41) or solving the generalized FKE (45) with a reflecting BC at the origin.One can appreciate that the latter task is more amenable.Making use of the method of images [60] the solution simply follows through knowledge of the Green's function of Eq. ( 45), that is the solution of the FDE with a reflecting BC, which is given in the Laplace domain by [1] W P α (x, λ α , t|x 0 ), is then constructed in terms of the Green's function as, From the defect technique [61] Eq. ( 48) can be written in the Laplace domain as follows [24], Substitution of Eq. ( 47) into Eq.( 49) gives, which can be readily shown to satisfy the radiation BC at x = 0.

First-Reaction Time Probability Density
To find the FRT probability density we first consider the survival probability Eq. ( 44), i.e.Q α (t|x 0 ) = P[T > t] x0 , to which we may write in terms of the survival probability of normal diffusion, As the FRT probability density, F α (t|x 0 ), is the density of the reaction times T , it is then related to the survival probability via F α (t|x 0 ) = −∂ t Q α (t|x 0 ).After using Eq. ( 11) we find which can be integrated by parts leading to where we have used g(τ, t) = δ(τ −S α (t)) .Comparing Eq. ( 53) to Eq. ( 19) one can see that the FRT probability density is simply related to This relation is obvious from the radiation BC since the FRT probability density is equal to the (negative) flux at the boundary.Then from Eq. ( 50) we find the Laplace Transform of the FRT density as [56], Since Eq. ( 55) is hard to invert directly, we proceed by finding the Mellin transform, M{f (t)} = f (s) = ∞ 0 t s−1 f (t)dt, and then performing the inverse transform [63,64] (see Appendix B and C).In the end we obtain, where Γ(z) = ∞ 0 t z−1 e −t dt is the Gamma function [65], γ(z, a) = a 0 t z−1 e −t dt is the lower incomplete Gamma function [65] and E a,b (z) is the two-parameter Mittag-Leffler function [66] (see Eq. (C11)).Using the integral form of the incomplete Gamma function, we may write Eq. ( 56) in a more compact integral form, (57) where φ(a, b; z) is the Wright function as defined by the series [67,68] for a > −1.From the Eqs.( 56) and ( 57) we can see that if the particle initially starts at the origin we obtain the FRT density in the simple form, Using the following Laplace transform relationship between the Wright function and the Mittag-Leffler function [69,70], for a > 0, then Eq. ( 57) becomes, This form of the FRT in Eq. ( 61) is useful for verifying the perfectly reacting/absorbing case (λ α → ∞) and the normal diffusion case (α → 1) (see Appendix D and E).Eq. ( 61) is now used to find the short and long time asymptotic form of the FRT.After a change of variable in the integral in Eq. ( 61) we use the expression to find the large t dependence since lower limit of the integral tends to zero.Thus at long times, using Eq. ( 60), we have the approximate form of F α (t|x 0 ), Using the asymptotic form of the Mittag-Leffler function [66] we have the following asymptotic dependence for t → ∞, which confirms the t −1−α/2 dependence [56] as well as that it possesses the same asymptotic dependence as in the perfectly reacting case [71], leading to an infinite mean [72].
For small t we the can disregard the integral on the right-hand side (RHS) of Eq. ( 66), then using the asymptotic form of φ(a, b, −z) for z → ∞ [68,73], we find F α (t|x 0 ) has the following dependence for t → 0, Eq. ( 67) indicates a short time exponential form for the FRT consistent with the perfectly reacting case [71].We plot F α (t|x 0 ) in Fig. (2) to prove the validity of our analytic results.

Survival Probability and Subordinated Local Time Density
Due to the reaction times, T , being governed by the subordinated local time, (S α (t)), we can find the distribution of (S α (t)), ρ α (l, t|x 0 ), and compare to
simulations.From Eq. ( 44) we know that ρ α (l, t|x 0 ) is the inverse Laplace transform of the survival probability Q α (t|x 0 ) with respect to λ α .By finding the survival probability we thus have an easy route to determine the subordinated local time density.We find Q α (t|x 0 ) = ∞ t F α (u|x 0 )du by using Eq. ( 61) and Eq. ( 65) along with the fact that [66] azφ(a, a to realize that From Eq. ( 69) we recover the known solutions for the perfectly reactive and normal diffusion cases (see Appendix D and E).Now we perfrom a change of variable on Eq. ( 69), which gives after integrating by parts and making use of Eq. ( 65) this then leads to From Eq. ( 71) the inverse Laplace transform is straightforward and we find the subordinated local time density to be, Clearly for a particle starting at the origin we obtain a simpler expression due to the first term on the RHS of Eq. ( 72) vanishing.
We are able to find the moments of ρ α (l, t|x 0 ), 44), where Q α (t|x 0 ) can be expressed as the moment generating function of (S α (t)), such that By repeatedly integrating by parts Eq. ( 69) whilst using Eq. ( 65), we are able to express Q α (t|x 0 ) as the following infinite series, Thus we find the moments are uniquely given in terms of Wright functions, For x 0 = 0 we obtain a simple form for the moments, which can be realised by using the series form of the Mittag-Leffler function in Eq. ( 59).Interestingly, for the specific value of α = 2/3 [69], we find the subordinated local time density to be in terms of the Airy function, Ai(z) = 1

Analytic Simulation
Fig. 3 The subordinated local time density, ρα(l, t|x 0 ), Eq. ( 77) for α = 2/3, plotted for different values of t, with x 0 = 0 and Kα = 1, with all quantities in arbitrary units.We plot these lines against stochastic simulations (shown as circles) of the subordinated local time, as described in Sec.4.4.
zt)dt for (z) = 0 [74], For α = 1, one can see that Eq. ( 72) reduces to the well known Gaussian solution (see Appendix E).In Figure (3), Eq. ( 77) is utilised to show the correct method of simulating the subordinated local time, as described in the next section.

Simulations
Due to the complexity of our analytic results it is important to validate with numerical simulations.A common method for simulating subdiffusion is through the CTRW formalism in the diffusive limit with a heavy-tailed waiting time distribution, however by keeping in the theme of this work we generate our simulations using a subordinated approach [45], as described in Fig. (1).Alternatively, one can simulate the subdiffusion trajectories via the algorithm in Ref. [75], which differs from Ref. [45] through not needing to explicity calculate S α (t).The aforementioned subordinated approach lends naturally to the simulation of the subordinated local time.We are able to effectively approximate, (S α (t)) via the discrete construction in Sec.2.3, where we approximate the limit in Eq. ( 25) for ε << 1.We approximate the integral in Eq. ( 24) by writing it in Riemann-Stieltjes form and using the left sided Riemann summation approximation, i.e. by discretizing the time interval, [0, t], as 0 = t 0 < t 1 < ... < t M = t, then As we are considering a bounded domain, the subdiffusive trajectories are made to be reflected at the origin.As pointed out in Ref. [20] one must ensure ε >> √ 2K α ∆t, where ∆t = t/M , such that the characteristic size of the jump length is not larger than the region ∂Ω, but the smaller the value of ε the better the approximation of the subordinated local time.The simulation for the FRT density naturally follows from the subordinated local time, where for each trajectory we draw a random number from an exponential distribution with mean 1/λ α and record the FRT when the subordinated local time exceeds the value of the random number.Given the good match between simulations and theory in Figs.(2) and ( 3) this provides validation for the analytic results presented.

Summary and Conclusions
In summary we have derived, through subordination techniques, a type of generalized FKE and used such an equation to understand and analyse subdiffusion in the presence of a reactive boundary.After deriving the α-stable inverse subordinator we perform a time change on the classical FKE to construct a generalized analogue, and find its solution.We introduce the notion of the subordinated local time and we interpret it as the continuum limit of the number of times a subdiffusive particle reaches a given location via a CTRW formulation.An important finding that emerges is that the time spent for each visit does not play a role in the subordinated local time.We apply the formalism to demonstrate how the generalized FKE can be used to describe a subdiffusion process in the presence of reactions.For that we show that the generalized FKE can be thought of as the position density of a subdiffusing particle with the requirement that the subordinated local time is less than a random variable drawn from an exponential distribution with a mean related to the parameter p in the generalized FKE.
We also consider what would be the relevant BC if the reaction occurred on a boundary.We study this aspect by considering a CTRW on a discrete lattice and introduce a reflection probability at the origin, k, which gives a general flux condition.From this condition we take the relevant limits to obtain the generalized form of the radiation BC associated with subdiffusion and find the generalized reactivity parameter, λ α .We then demonstrate the equivalence between the generalized FKE with a reflecting boundary and the FDE with a radiation boundary.
We employ the generalized FKE to study the FRT of a subdiffusive particle in the presence of a radiation boundary at the origin.We solve the generalized FKE using Green's function techniques and obtain an analytic solution in the Laplace domain.We then find the relation between the solution of the generalized FKE and the FRT probability density and obtain this quantity in the Laplace domain, which we are able to invert by converting it into a Mellin transform.The FRT probability density is obtained in terms of the Mittag-Leffler function and an infinite series of lower incomplete Gamma functions or alternatively as an integral involving the Wright function.From this we are able to analyse the short and long time asymptotic form of this density recovering expected dependencies.Due to the fundamental connection between the FRT and subordinated local time, we calculate the subordinated local time density and all its moments.Finally, we show how our analytic results match with simulations, proving the validity of simulating subdiffusion in the presence of a radiation boundary using the subordinated local time approach.
A natural extension to this work would be to consider the subordinated occupation time functional (U (x) = I Ω (x), for some spatial region Ω) and how this may be used as in a similar sense as the subordinated local time here to describe subdiffusion with reactions in a certain region, not just at a boundary.Further future directions could include developing the backward version of the generalized FKE considered here to study the density of various other subordinated functionals and look at how they compare to functionals of subdiffusion [9].Clearly functionals dependent on the underlying subdiffusive path like local time, occupation time etc. are certainly going to be different.However, the generalized FKE may be applicable to a class of functionals associated with first-passage times, due to only needing the knowledge of whether the particle has reached a specific point, not how long it has been there.leading to From the series representation of the Mittag-Leffler function [66],  then after writing back in terms of γ(z, a) we recover the summation in Eq. (56).Note this sum should start from n = 1 due 1/Γ(0) = 0 for n = 0, but we leave it in this form for easier relation to the Wright function.