# Stochastic modeling of movement for Helium particles in a graphite channel

Jingqi Zhang, Mitra Fouladirad, Nikolaos Limnios, Pierre Magnico

► To cite this version:

Jingqi Zhang, Mitra Fouladirad, Nikolaos Limnios, Pierre Magnico. Stochastic modeling of movement for Helium particles in a graphite channel. Physica A: Statistical Mechanics and its Applications, 2025, 675, pp.130818. 10.1016/j.physa.2025.130818. hal-05246509

HAL Id: hal-05246509

https://hal.science/hal-05246509v1

Submitted on 9 Sep 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Stochastic modeling of movement for Helium particles in a graphite channel

Jingqi Zhang<sup>a,c,*</sup>, Mitra Fouladirad<sup>b,●</sup>, Nikolaos Limnios<sup>c</sup>, Pierre Magnico<sup>b</sup>

<sup>a</sup> LIST3N, Université de Technologie de Troyes, Troyes, 10004, France

<sup>b</sup> M2P2, Aix Marseille Université, CNRS 7340, Centrale Med, Marseille, 13013, France

<sup>c</sup> LMAC, Université de Technologie de Compiègne, Compiègne, 60319, France

## ARTICLE INFO

### Keywords:

- Knudsen diffusion

- Semi-Markov process

- Asymptotic result

- Diffusion coefficient

## ABSTRACT

In this article, we present a stochastic model for the movement of Helium particles within a graphite channel, focusing on Knudsen diffusion. We develop a semi-Markov model to describe the movement of the particle, derive the stationary distribution of its mean position, and analyze the model's asymptotic properties. To validate the model, we compare its theoretical outcomes with Monte Carlo simulations. As temperature significantly influences on the movement of particles, two situations are studied for high and low temperature. In both cases, theoretical and simulation results by Monte Carlo coincide. Furthermore, we propose estimation methods for the local parameters of the model and demonstrate its application using data from Molecular Dynamics simulations.

## 1. Introduction

Rarefied gas flows are encountered in numerous industrial fields such as gas storage [1–3], micro-vacuum technologies [4], advanced material synthesizing [5,6], membrane separation [7,8], fuel cell [9,10], gas sensors [11,12], shale gas recovery [13,14] and also cement-based and concrete materials deterioration [15,16]. In nanochannels, the gas mean free path length $l$ is of the order or greater than the channel width $H$ and gas/surface collisions govern the mass transfer. At this scale, the continuum assumption of flows is no more valid and in the free molecule regime, i.e., when $K_n = l/H \gg 1$, the molecule trajectory is a series of free-flights interrupted by collisions with channel walls. Therefore, it is of crucial importance to understand and predict the gas/wall momentum and the energy transfer and also to quantify their influence on mass transfer in nanometric systems.

The gas/wall boundary condition is generally described by a scattering kernel which gives a statistical relation between the incident velocity and the reflected one after the gas molecule strikes the solid surface and leaves the near-wall layer. The Maxwell [17] and Cercignani–Lampis [18,19] models are mainly used owing to the simplicity of their formulation. Indeed, their scattering kernel is governed by one or two empirical parameters named accommodation coefficients (AC) which describe how the kinematic properties of the incident gas molecules accommodate to the surface. It also reveals the correlation between the incident and the outgoing velocity. The main drawbacks of these models are the independence of AC from the molecular velocity components and from the incidence angle, and the lack of consideration of the gas/wall physico-chemical interactions which control the surface diffusion.

These interactions are characterized by a potential well which depth $W$ results not only in the contribution of all the solid atoms but also in the contribution of their local configuration at the surface, i.e., the surface roughness [20,21]. Depending on its kinetic
energy, a gas molecule can bounce several times on the wall before leaving the surface when the normal kinetic energy becomes greater than the local potential energy. Therefore, the molecule can spend a time and displace along the surface. This mechanism, called surface diffusion [22–24], must be taken into account when the width of the near wall resident time distribution (DTS) is greater than the mean free-flight time between two collisions.

To improve the interplay between the scattering process and the microscopic gas/wall interaction, several physical models have been proposed such as hard/soft cube model, washboard model or Yamamoto model among many others [25–27,27,28]. Some authors also used with machine learning statistical approach to determine a scattering kernel in order to couple better the collision process with the incident condition or to determine a fully statistical kernel [29,30]. As concerns the washboard model, improvements have been made to solve the problem of reciprocity violation and to make the model free of unknown parameters in the case of light gas atoms colliding heavy solid atoms [31,32]. These improvements make the model promising despite the complexity of the final kernel expression.

In connection to the present work, Knudsen self-diffusion in a 3D slit pore was first studied by molecular dynamics (MD) in [33]. To simplify the diffusion process, the He/graphene system is chosen, the gas/gas interaction is neglected, i.e the Knudsen number $Kn$ is infinite, and the system is in thermal equilibrium. To link the self-diffusion to the gas/wall collision process at an imposed temperature, the tangential momentum accommodation coefficient (TMAC) and the near wall DTS were determined from molecular dynamics computation consisting in throwing He atoms to the wall from a transversal distance $z_{cut}$ which is the He/C 12-6 Lennard Jones potential cutoff distance and characterizes the depth of the near wall region. For the self-diffusion investigation, the two longitudinal velocity components and the transversal one follow the Maxwellian and the Rayleigh one respectively with a null mean value owing to the thermal equilibrium. As concerns the collisional investigation, these velocity distributions are imposed at the inlet of the near wall region, but the Rayleigh distribution of the transversal component has a non null mean value because the gas atoms are thrown to the wall. The two investigations are carried out with thermostatic walls composed of three graphene sheets. Molecular dynamics shows that at high temperature (T > 200K), the collision process is mostly specular, i.e. the TMAC is far smaller than one. The near wall resident time distribution is a combination of a Fréchet distribution of order 2 at short time, in accordance with the Rayleigh distribution of the normal gas velocity, and an exponential one at long time which time duration and time decrement decrease with the increasing Temperature. But at small temperatures (T < 100K), the collision process is governed by the surface diffusion which induces a high accommodation and a diffusive gas reflection with a TMAC value close to one. In addition, a deviation with order 2 was observed and can be arisen from the fact that the total energy of the gas atoms becomes comparable to the potential energy, and as the gas atoms approach the wall, their velocity increases significantly. In the case of hard spheres diffusing in the rough slit pore, authors extended the Smoluchowski model by using the CL model instead of the Maxwell one [34,35]. This model agrees with the molecular dynamics results at high temperature but fails at small temperature. Under low temperature, a Monte Carlo approach including the order 2 Fréchet distribution associated with an exponential distribution, the TMAC, and the surface diffusion is proposed [33]. However, it is necessary to improve the model by using generalized Fréchet distribution to better present DTS at low temperature.

To further refine the understanding of Knudsen diffusion in nanometric systems and address the limitations of existing models at low temperatures, the adoption of a stochastic framework emerges as a natural and effective solution. Stochastic modeling is particularly well-suited to capturing the inherent randomness in gas/surface interactions, as it allows for the inclusion of probabilistic behaviors, such as the near-wall DTS and its temperature dependence. We adopt a continuous-time framework, as it aligns naturally with the real-time evolution of particle dynamics, where events and transitions occur continuously rather than in discrete intervals. Furthermore, a countable discrete state space is employed, balancing mathematical simplicity with the capacity to accurately describe the system. A non-countable state space would introduce unnecessary complexity without offering additional insights for the current scope of analysis. The Markov property, which assumes that the future evolution of the system depends solely on its current state, offers significant advantages for modeling gas dynamics, where memoryless state transitions are valid. While Markov processes provide a powerful modeling approach, their assumption of exponentially distributed residence times limits their applicability in systems where transitions depend on more general distributions.

The semi-Markov Process (SMP) has property that the sojourn time of the current state depends only on the present state, without being constrained by exponential distribution [36]. It generalizes Markov and renewal processes by allowing the residence time distribution to take any form on $R_+$. In [37], the authors compare the evaluation of machinery readiness using a Markov process and a SMP, finding that the Markov model has a result 46% lower than that of the SMP. Therefore, in situations where the sojourn time is non-exponential, the SMP offers greater flexibility and more accurate results. As in our case, the sojourn time follows a mixture of a Fréchet distribution and an exponential distribution, SMP is more suitable. Besides, this flexibility has led to its successful application across diverse fields, including economic forecasting [38], epidemic modeling (e.g., SARS) [39], earthquake prediction [40], queuing theory [41], DNA analysis [42], reliability analysis [43], and open quantum systems [44].

In this study, we introduce for the first time a semi-Markov process model to analyze Knudsen diffusion in a 2D slit pore system, focusing on the critical moments when helium particles enter and exit the interaction zone. This approach builds on the stochastic relationship introduced in [33], combining analytical asymptotic results with numerical examples and comparisons to Monte Carlo (MC) simulations. Through this model, we derive key theoretical results, including the stationary distribution for the mean residence time of the system, the velocity distribution function in the x-direction, and asymptotic expressions for the time distribution and the x-direction travel distance. Furthermore, a closed-form expression for the diffusion coefficient is obtained, significantly reducing the computational time compared to direct simulation methods, and establishing a relationship between the diffusion coefficient and the model parameters. As previously highlighted, the DTS transitions from a Fréchet distribution for short residence times to an exponential decay at longer times. The transition between these behaviors, marked by a time threshold $t_0$, is challenging to identify.
To address this, we propose several methods for time threshold and parameter estimation, which are compared with MD results. Nevertheless, the current models do not fully capture the collision mechanism, as the primary mode of mass transfer is governed by diffusion, which plays a dominant role in this system.

The remainder of the paper is organized as follows. In Section 2, the physical model is detailed. The expression of the resident time close to the wall, the velocity and the position of particles are given. In Section 3, we propose an extended semi-Markov process to model the particle behavior. The stationary distribution is derived and the asymptotic behavior is studied. In Section 4, we additionally carry out MC simulations although more time consuming it permits to compare and validate the analytical results. We use statistical methods to estimate time threshold and parameters. The methods are first validated using MC simulation data and subsequently applied to MD data to derive the theoretical results. Section 5 provides conclusions and outlines perspectives for future research. Finally, the proofs of theoretical results, along with additional details on optional estimation methods and their verification, are presented in the Appendices.

### Notation.

$m_0$ particle mass.

$d$ parallel wall distance.

$z_{cut}$ absorption distance for interaction.

$v_z$ velocity in z-direction.

$v_x$ velocity in x-direction after collision at wall.

$v'_x$ velocity in x-direction before collision at wall.

$v_x^M$ Maxwellian distribution.

$T$ temperature.

$\alpha_{tan}$ tangential accommodation coefficient.

$k_B$ Boltzmann constant.

$t_0$ time threshold of resident time at wall between long resident regime and short resident time regime.

$D_s$ diffusion coefficient of surface.

$D_{channel}$ diffusion coefficient of channel.

$J_n$ Embedded Markov Chain for jump time.

$f_N$ density distribution function of resident time at wall.

$F_N$ cumulative distribution function of resident time at wall.

$V_z^{(n)}$ velocity in z-direction after n-th change of state, r.v.

$V_x^{(n)}$ velocity in x-direction after n-th change of state, r.v.

$N_t$ number of change of state before time t.

$S_n$ time of the n-th change of state.

$X_n$ resident time between (n-1)th and n-th change of state.

$a$ travel distance function between (n-1)th and n-th change of state.

$Q$ semi-Markov kernel.

$P_t$ transition matrix function of semi-Markov process.

$\epsilon_z$ minimum velocity via z-direction.

$\bar{a}_s$ convergence almost surely.

$\bar{P}$ convergence in probability.

$\bar{D}$ convergence in distribution.

We denote by uppercase letters the random variable (r.v.) and by lowercase letters their values, i.e., $V_x$ is a r.v. and $v_x$ its value.

## 2. Physical model

In this section, we provide a detailed description of the physical system, presented in [33], under consideration, which serves as the basis for the stochastic modeling approach introduced in the following section. This model is divided into several distinct phases, each of which is carefully analyzed in the subsections that follow. First, the particle's velocity distribution is considered, with particular focus on the velocity components in the longitudinal and transversal directions. The travel distance along the longitudinal direction is presented, followed by a detailed examination of the residence time distribution, which plays a crucial role in the stochastic description of the system's dynamics. By integrating these components, the model provides a comprehensive framework for understanding the transport properties of helium gas in nanoscale systems and serves as a foundation for the subsequent mathematical modeling.

### 2.1. Model description

The model involves a system in thermal equilibrium with a monoatomic gas, specifically Helium. The Helium gas moves freely with a null mean longitudinal velocity within a graphite channel characterized by a larger width, denoted by $d$. The channel's
<figure><img src="image_2.png" /><figcaption>Fig. 1. System collision process.</figcaption></figure>

surface is notably smooth. In this system, interactions between gas/gas are considered negligible due to the nanoscale dimension of the channel, resulting in a gas mean free path that exceeds the channel size. Consequently, we can focus on the trajectory of individual Helium particles within the channel. Given the scarcity of particles moving perfectly to the channel, we assume that each particle collides with the channel with an angle, while the channel itself remains immobile. The collision process is illustrated in Fig. 1 and can be divided into two distinct phases: the free-flight phase and the interaction phase. During the free-flight phase, the particle's velocity is assumed to remain constant. The resident time of this phase is determined by the velocity in z-direction, denoted by $v_z$. During this phase, the particle has a free-flight distance, labeled as $h$, with a mass of $m_0$. When the particle reaches a specific distance from the surface, defined as $z_{cut}$, the interaction phase is initiated. In interaction phase, the particle is absorbed by the surface, stays for a resident time $\tau$, and is eventually desorbed to proceed into another free-flight phase. As the two channel walls are parallel, we have specifically $h = d - 2z_{cut}$.

We focus on a particle within this system which is simplified in 2D as the walls are isotropic in the longitudinal direction and for this reason the particle behaves the same way in direction x and y. We suppose that the particle possesses an initial position $(x_0, z_0)$ and an initial velocity $V^{(0)} = (V_x^{(0)}, V_z^{(0)})$. It passes through a channel with distance $d$. The velocity of the nth change of zone is denoted by $V^{(n)} = (V_x^{(n)}, V_z^{(n)})$. The distance realized toward the x axis of collision is denoted by $\Delta x$. We use the velocity vector $V = (V_x, V_z)$.

Under high temperature, the reflection is quasi instantaneous. Let us denote by $V_x^{pre}$ the pre-collision velocity and $V_x^{post}$ the post-collision velocity. For each collision, the velocities exhibit a relationship described by the following formula [33]:

$$V_x^{post} = (1 - \alpha_{\text{tan}})V_x^{pre} + \sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}V_x^M \quad (1)$$

Here, $\alpha_{\text{tan}}$ is the TMAC, which depends on temperature $T$. TMAC is related to the fraction of diffuse collisions, where the reflected velocity is uncorrelated with the incident velocity. By neglecting the gas/wall Lennard-Jones interaction, the local variation of the surface orientation increases with the roughness amplitude, causing TMAC to approach 1. This indicates that the reflected velocity becomes increasingly decorrelated from the incident velocity. In the context of graphene, the surface roughness is extremely low, which implies that at high temperatures, TMAC remains close to 0 (approximately 0.04 for helium), which means that the collisions are nearly specular. Various models relating TMAC to surface roughness have been referenced in the extended introduction of the previous version. $V_x^M$ follows the Maxwellian distribution, represented by a Gaussian distribution $V_x^M \sim N(0, k_B T/m_0)$, where $k_B$ is the Boltzmann constant.

In this context, it is relevant to consider that surface velocity of atoms remains constant. The surface velocity is approximately equal to:

$$V_{sx} = \frac{1}{2}(V_x^{post} + V_x^{pre}).$$

Therefore, the ballistic displacement is given by:

$$(\Delta x)_{surf} = V_{sx} \tau. \quad (2)$$

Under low temperature, the behavior of the system varies based on the resident time. For short resident times ($\tau < t_0$), the collision process exhibits ballistic reflection characteristics, following the same relation as Eq. (2). However, for longer resident time ($\tau \ge t_0$), the collision behavior comprises two stages: the ballistic stage and the diffusive stage, i.e.,

$$(\Delta x)_{surf} = (\Delta x)_{bal} + (\Delta x)_{dif}.$$

During the ballistic stage before $t_0$, the displacement, represented as $(\Delta x)_{bal}$, and we have,

$$(\Delta x)_{bal} = t_0 V_{sx}, \quad V_{sx} = \frac{1}{2}(V_x^{post} + V_x^M).$$
In the diffusive stage, that occurs after $t_0$, the displacement $(\Delta x)_{dif}$ follows a Gaussian distribution given by:

$$ (\Delta x)_{dif} \sim \mathcal{N}(0, 2D_s(t - t_0)), $$

where $D_s$ is the diffusion coefficient, which depends on the temperature. This coefficient $D_s$ represents the rebounds of diffusion on surface. The rebounds occur in random directions. After a sufficiently large number of rebounds, the surface transport process becomes diffusive, and the surface mean squared distance (MSD) increases proportionally to $t$ when the correlation time is sufficiently large.

The velocity component in the z-direction $V_z$, follows a Rayleigh distribution, and this characteristic remains consistent across all temperature conditions.

$$ V_z \sim \mathcal{R}(0, \theta^2), $$

where $\theta = \sqrt{k_B T / m_0}$ is the reference thermal velocity.

## 2.2. Resident time at wall layer

In [33], the MD shows that the transport of Helium particle is essentially ballistic under high temperature ($T > 200K$), and two modes of transport, i.e., ballistic and diffusive, exist with long resident time under low temperature ($T < 100K$). With short resident time, the model in [33] suppose that Helium particles leave the interaction zone with a single bounce. As the velocity in z-direction follows Rayleigh distribution, the resident time is proportional with the inverse of $V_z$ for short resident time.

Here, we assume that the resident time at wall layer follows Fréchet distribution $F_r(s, \beta, m)$, where $s$ is the shape parameter, $\beta$ the scale parameter, and $m$ the location of minimum. Let $f_0(t)$ denote the probability density function (pdf) of Fréchet distribution $F_r(s, \beta, m)$, i.e.,

$$ f_0(t) = \begin{cases} 0, & t < m, \\ \frac{s}{\beta} \left(\frac{t-m}{\beta}\right)^{-1-s} e^{-\left(\frac{t-m}{\beta}\right)^{-s}}, & m \le t, \end{cases} $$

and $F_0(t)$ its cumulative distribution function (cdf),

$$ F_0(t) = \begin{cases} 0, & t < m, \\ e^{-(t/m)^{-s}}, & m \le t. \end{cases} $$

When it comes to long resident time, the model suppose that Helium particles leave the interaction zone with several bounces, and exchange of energy and movement quantity, which result in the resident time following an exponential law. Let us denote by $f_1(t)$ the shifted exponential pdf, i.e.,

$$ f_1(t) = \lambda e^{-\lambda(t-t_0)}, \quad t \ge t_0, $$

with $\lambda > 0$. And $F_1(t)$ the corresponding cdf,

$$ F_1(t) = 1 - e^{-\lambda(t-t_0)}. $$

The pdf $f_T(t)$ of resident time is expressed using the weight $p$, which is the probability of entering short resident regime, and $q = 1 - p$ the probability of entering long resident regime, given by:

$$ f_T(t) = \begin{cases} \frac{p}{f_0(t_0)} f_0(t), & t < t_0, \\ q f_1(t - t_0), & t_0 \le t, \end{cases} \\ = \begin{cases} 0, & t < m, \\ \frac{p}{e^{-(t_0/m)^{-s}} \beta} \frac{s}{\beta} \left(\frac{t-m}{\beta}\right)^{-1-s} e^{-\left(\frac{t-m}{\beta}\right)^{-s}}, & m \le t < t_0, \\ q \lambda e^{-\lambda(t-t_0)}, & t_0 \le t \end{cases} \tag{3} $$

with $m \ge 0$, and $p$ the weight parameter, characterized by model parameters and fixed transition time $t_0$. Therefore, the corresponding cdf $F_T$ is normalized and given by:

$$ F_T(t) := \int_0^t f_T(x) dx = \begin{cases} 0, & t < m, \\ \frac{p}{e^{-(t_0/m)^{-s}} \beta} e^{-\left(\frac{t-m}{\beta}\right)^{-s}}, & m \le t < t_0, \\ p + q(1 - e^{-\lambda(t-t_0)}), & t_0 \le t. \end{cases} \tag{4} $$

Here, the pdf $f_T$ is not continuous, while the cdf $F_T$ is continuous but not differentiable at $t_0$.
## 2.3. Velocity

The velocity in z-direction, after the nth collision, $V_z^{(n)}$ for $n \in \mathbb{N}$, are identically and independently distributed (i.i.d) which follow Rayleigh distribution function $\mathcal{R}(0, \theta^2)$, with pdf $f_R$ as follows,

$$f_R(v_z) = \frac{v_z}{\theta^2} e^{-\frac{v_z^2}{2\theta^2}},$$

where $\theta = \sqrt{k_B T / m_0}$.

Let us denote by $f_N$, the probability density function associated to a centered Gaussian random variable with variance $\theta^2$, $\mathcal{N}(0, \theta^2)$, as follows:

$$f_N(v) = \frac{1}{\theta \sqrt{2\pi}} e^{-\frac{1}{2} \frac{v^2}{\theta^2}}.$$

The velocity in x-direction depends on the temperature, the resident time and its previous velocity as described in Eq. (1) for $n \in \{1, 2, ...\}$. To account for these factors, we consider the transition kernel of the velocity under:

* High temperature

$$P(u, dv) = \mathbb{P}(V_x^{(n)} \in dv | V_x^{(n-1)} = u) = f_N\left(\frac{v - (1 - \alpha_{\text{tan}})u}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}},$$

* Low temperature

$$P(u, dv) = \begin{cases} \mathbb{P}(V_x^{(n)} \in dv | X_n < t_0, V_x^{(n-1)} = u), \\ \mathbb{P}(V_x^{(n)} \in dv | X_n \ge t_0, V_x^{(n-1)} = u) \end{cases} \\ = \begin{cases} f_N\left(\frac{v - (1 - \alpha_{\text{tan}})u}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}, \\ f_N(v)dv. \end{cases}$$

with $V_x^{(0)} \sim \mathcal{N}(0, \theta^2)$ and $X_n$ the corresponding resident time at wall layer.

## 2.4. Travel distance

The distance realized by the particle depends on temperature, velocity and resident time at wall. We also consider the distance under two situations:

* High temperature

$$\Delta x = \frac{1}{2}(v_x^{\text{pre}} + v_x^{\text{post}})t,$$

* Low temperature

$$\Delta x = \begin{cases} \frac{1}{2}(v_x^{\text{pre}} + v_x^{\text{post}})t, & t < t_0, \\ \frac{1}{2}(v_x^{\text{pre}} + v_x^{\text{post}})t_0 + y, & t \ge t_0, \end{cases}$$

with y the value, and its r.v., $Y \sim \mathcal{N}(0, 2D_s(t - t_0))$ when $t \ge t_0$.

## 3. Stochastic modeling

In Section 2.2, we introduced that the DTS follows exponential distribution for long time regime, while it follows Fréchet distribution for short time regime. This complexity necessitates a more flexible framework than a Markov process. The SMP provides such an extension by relaxing the constraint of exponentially distributed waiting times between transitions, allowing for arbitrary residence time distributions. This makes the SMP well-suited for capturing deviations observed in the DTS. By employing this approach, it becomes possible to reconcile the MD results with theoretical predictions and to derive generalized distributions for the DTS as well as refined estimates for diffusion coefficients. These theoretical insights are critical for understanding the interplay between Knudsen and surface diffusion, enabling more accurate predictions of gas transport properties in nanometric systems.

In this section, we develop the stochastic modeling framework described in Section 2, utilizing the SMP. The SMP is defined on a discrete state space and in continuous time, providing a flexible framework for modeling the dynamics of the whole system, including free-flight zone and interaction zone. In this system, we define state 1 the free-flight phase and state 2 the interaction phase. In state 1, the velocity is constant, and the resident time in state 1 depends only on the velocity in z-direction, $h/V_z^{(n-1)}$. State 2 is the collision of particle, where it hits on wall and changes its velocity. The resident time follows the distribution $F_T$. From this SMP model, we derive the stationary law governing the mean residence time and the velocity distribution in the x-direction. The Law of Large Numbers (LLN) and central limit theorem (CLT) are then applied to show the convergence of the time distribution and the mean travel distance along the x-direction [45]. Finally, we derive the diffusion coefficient, providing key theoretical results for the system.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Resident time representation.</figcaption></figure>

## 3.1. Semi-Markov process

According to semi-Markov theory [36], defining the SMP and its transition matrix are crucial for obtaining the asymptotic results. We now define a semi-Markov process ($Z_t$)$_{t \ge 0}$ with state space $E = \{1, 2\}$, and the embedded Markov Renewal Process (MRP) ($J_n$, $S_n$)$_{n \in \mathbb{N}}$, where ($J_n$)$_{n \in \mathbb{N}}$ is the Embedded Markov Chain (EMC), and $S_n$ the nth jump of the SMP, $0 = S_0 < S_1 < \cdots < S_n < \cdots$. The resident time in each state $X_n = S_n - S_{n-1}$. Their relation is presented in Fig. 2.

The related counting process ($N_t$)$_{t \in \mathbb{R}_+$} [36] is defined as

$$N_t = \sup\{n \ge 0 : S_n \le t\}.$$

The resident time in state 1 follows the distribution function $G(t)$,

$$G(t) = \mathbb{P}\left(\frac{h}{V_z^{(n-1)}} \le t\right) = e^{-\frac{h^2}{2\theta^2}} \quad t > 0.$$

The resident time in state 2 follows the distribution function $F_T$, see Eq. (4). Therefore, the semi-Markov kernel is given by:

$$Q(t) = \begin{pmatrix} 0 & G(t) \\ F_T(t) & 0 \end{pmatrix},$$

which represents the probability that during a time interval $(0, t]$ the system jumps from state $i$ to state $j$, $i, j \in E = \{1, 2\}$.

As it is mentioned by [36], the transition matrix $P_t$ associated to the kernel $Q(t)$ is then given by:

$$P_t = M * \begin{pmatrix} \bar{G} & G * \bar{F}_T \\ \bar{G} * F_T & \bar{F}_T \end{pmatrix},$$

where $M = 1 + \sum_{n \ge 1} (F_T * G)^{(n)}$, $\bar{G} = 1 - G$, $\bar{F}_T = 1 - F_T$.

The symbol * here stands for the Stieltjes convolution on $\mathbb{R}_+$ defined by:

$$F_T * G = \int_0^t F_T(t-s)dG(s),$$

and the power $F_T^{(n)}$ means the n-fold convolution of $F_T$, i.e., $F_T^{(n)} = F_T * \cdots * F_T$, n times.

The Embedded Markov Chain has the transition matrix defined by $\mathbf{P} = \lim_{t \to \infty} Q(t)$, which gives

$$\mathbf{P} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix},$$

which signifies that the system surely transits from free-flight phase to interaction phase or inversely.

## 3.2. Travel distance

Our objective is to investigate the properties of the travel distance and subsequently deduce the diffusion coefficient of channel, so we define a function $a$ in this extended state space: $E \times \mathbb{R}_+ \times \mathbb{R} \times \mathbb{R} \to \mathbb{R}$ for the travel distance in each state.

(a) Under high temperature

$$a(i, t, u, v) \equiv a_i(t, u, v) = \begin{cases} ut, & i = 1, \\ \frac{u+v}{2}t, & i = 2. \end{cases}$$
(b) Under low temperature

$$a(i, t, u, v) = \begin{cases} ut, & i = 1, \\ \frac{u+v}{2}t, & t \le t_0, \\ \frac{u+v}{2}t_0 + y, & t > t_0, \end{cases}$$

where y is the value of r.v. of $Y \sim \mathcal{N}(0, 2D_s(t - t_0))$.

Now we consider an embedded MRP ($J_n, X_n, V_x^{(n-1)}, V_x^{(n)}$)$_{n \ge 0}$ corresponding to travel distance.

The average travel distance until the nth change of state is defined by:

$$A_n := \sum_{i=1}^{n} a(J_i, X_i, V_x^{(i-1)}, V_x^{(i)})/n.$$

The average travel distance up to time t, is:

$$A_t := \left( \sum_{n=1}^{N_t} a(J_n, X_n, V_x^{(n-1)}, V_x^{(n)}) + a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)}) \right) / t.$$

### 3.3. Embedded Markov renewal process

Let us define also the indicator function,

$$\mathbb{1}_A(x) = \begin{cases} 1, & \text{if } x \in A, \\ 0, & \text{if } x \notin A. \end{cases}$$

When A has only one element {$v$}, we denote $\mathbb{1}_{\{v\}}(dv')$ as $\mathbb{1}_{v}(dv')$.

The transition matrix associated to the embedded MRP for the domain of travel distance is denoted as $\tilde{P}$, given by

$$\tilde{P}(i, t, u, v; di', dt', du', dv')$$

$$= \mathbb{P}(J_{n+1} = i', X_{n+1} \in dt', V_x^{(n)} \in du', V_x^{(n+1)} \in dv' | J_n = i, X_n = t, V_x^{(n-1)} = u, V_x^{(n)} = v)$$

$$= P_{ii'} F_i(dt') \mathbb{1}_v(uv') P_i(uv, dv') = F_i(dt') \mathbb{1}_v(uv') P_i(uv, dv').$$

Here we consider $P_i(u, dv)$ the transition of velocity given the state $i$, which depends on the temperature.

(a) Under high temperature

$$P_i(u, dv) = \begin{cases} \mathbb{1}_u(uv), & i = 1, \\ f_N\left(\frac{v(1-\alpha_{\text{tan}})u}{\sqrt{\alpha_{\text{tan}}(2-\alpha_{\text{tan}})}}, \frac{dv}{\sqrt{\alpha_{\text{tan}}(2-\alpha_{\text{tan}})}}\right), & i = 2. \end{cases}$$

(b) Under low temperature

$$P_i(u, dv) = \begin{cases} \mathbb{1}_u(uv), & i = 1, \\ f_N\left(\frac{v(1-\alpha_{\text{tan}})u}{\sqrt{\alpha_{\text{tan}}(2-\alpha_{\text{tan}})}}, \frac{dv}{\sqrt{\alpha_{\text{tan}}(2-\alpha_{\text{tan}})}}\right), & t \le t_0, \\ f_N(uv)dv, & t > t_0, \end{cases}$$

### 3.4. Stationary distribution

Here, we derive the stationary distribution governing the system, including the stationary distribution for the EMC, SMP, and the extended MRP. Additionally, we establish the stationary distribution of the particle's velocity in the x-direction, providing a comprehensive characterization of the system's long-term behavior.

The stationary law of the EMC derived from transition matrix P is:

$$v = \left(\frac{1}{2}, \frac{1}{2}\right). \qquad (5)$$

The stationary law $\pi$ of the SMP, which presents the probability that the system in state $i$, is stated as follows, see e.g. in [36]:

$$\pi_i = \frac{v_i m_i}{m_s},$$

where $m_s$ the mean resident time of the system.

In this case, the mean resident time in each state is given by:

$$\begin{align} m_1 &= \int_0^\infty \bar{G}(t)dt = \frac{h}{\theta} \sqrt{\frac{\pi}{2}}, \tag{6} \\ m_2 &= \int_0^\infty \bar{F}_T(t)dt = \frac{p}{e^{-(t_0/m)^{-s}} - 1} \int_0^{t_0} tdF_0(t) + q(t_0 + \frac{1}{\lambda}). \nonumber \end{align}$$
Therefore,

$$ \pi = \left( \frac{m_1}{2m_s}, \frac{m_2}{2m_s} \right), \qquad (7) $$

with

$$ m_s = \frac{1}{2}m_1 + \frac{1}{2}m_2. \qquad (8) $$

We define the stationary distribution $\tilde{\pi}$ of $(J_n, X_n, V_x^{(n-1)}, V_x^{(n)})$ [36]. Especially for simplification, we put

$$ F_{X_i}(t) = \begin{cases} G(t), & i = 1, \\ F_T(t), & i = 2. \end{cases} $$

Therefore, the stationary distribution of the extended MRP is,

$$ \tilde{\pi}_i(t, u, v) = \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} v_i F_{X_i}(dt) \rho(u) du P_i(u, dv), $$

where $\rho$ the stationary density of $V_x^{(n)}$.

Now, we try to find the stationary distribution for the velocity in x-direction. Let $F_\rho$ be the stationary distribution function on $\mathbb{R}$ which should satisfy

$$ F_{\rho}(B) = \int_{\mathbb{R}} F_{\rho}(du) \int_{B} f_{N}\left(\frac{v - (1 - \alpha_{\text{tan}})u}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}. \qquad (9) $$

**Proposition 1.** With initial x-direction velocity $V_x^{(0)} \sim \mathcal{N}(0, \theta^2)$, we have

$$ \mathbb{E}(V_x^{(n)}) = 0, \quad \text{Var}(V_x^{(n)}) = \theta^2. $$

**Proof.** See Appendix A. $\square$

**Theorem 1.** The stationary pdf $\rho$ of the Markov Chain $(V_x)$ is the unique solution of Eq. (9):

$$ \rho(x) = \frac{1}{\theta\sqrt{2\pi}} e^{-\frac{1}{2}\frac{x^2}{\theta^2}} \sim \mathcal{N}(0, \theta^2). $$

**Proof.** See Appendix B. $\square$

By analyzing the stationary distribution of the x-direction velocity, we gain insight into the long-term behavior of the particles, which reveals that the particles' x-direction velocities follow a Gaussian distribution in the long run.

### 3.5. Model asymptotic behavior

Here, we establish key theoretical results for the system's resident time and average travel distance using rigorous stochastic analysis. By using the Strong Law of Large Number (SLLN) [36] and the Central Limit Theorem (CLT) [46], we demonstrate the convergence properties of the resident time and derive asymptotic behavior for the particle's motion, ensuring a robust understanding of the system's long-term dynamics of the whole system, including free-flight phase and interaction phase. Eventually, we deduce the diffusion coefficient of channel.

In order to simplify the formulas, we use the following notations in the sequel:

* Exponential integral

$$ Ei(x) = - \int_{-\infty}^{\infty} \frac{e^{-t}}{t} dt = \int_{-\infty}^{x} \frac{e^t}{t} dt, \quad t \neq 0, \quad x < 0. $$

$$ E_1(x) = \int_x^\infty \frac{e^{-t}}{t} dt, \quad t \neq 0, \quad x > 0. $$

* Error function

$$ \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt. $$
These functions have corresponding table to consult for numerical values. Other terms not specified here are supposed to be calculated numerically.

It is well known from the Markov renewal theory that, see e.g., [39]:

$$ \frac{1}{t} N_t \underset{t \to \infty}{\stackrel{\text{a.s.}}{\longrightarrow}} \frac{1}{m_s}, \qquad (10) $$

where $m_s = (m_1 + m_2)/2$, which is stated in Eq. (8).

Moreover,

$$ \frac{N_t - t/m_s}{\sqrt{t\sigma^2/m_s^3}} \underset{t \to \infty}{\stackrel{D}{\longrightarrow}} \mathcal{N}(0, 1). $$

**Lemma 1.** We have the variance of the resident time

$$ \sigma^2 = \frac{1}{2}(\mathbb{E}(X_1^2) + \mathbb{E}(X_2^2)) - \frac{1}{4}(m_1 + m_2)^2, $$

with

$$ \begin{aligned} \mathbb{E}(X_1^2) &= -\frac{h^2}{2\theta^2} Ei\left(-\frac{\epsilon_z^2}{2\theta^2}\right), \\ \mathbb{E}(X_2^2) &= \frac{p}{e^{-(t_0-m)/\beta}} \int_0^{t_0} t^2 dF_0(t) + q(t_0^2 + \frac{2t_0}{\lambda} + \frac{2}{\lambda^2}), \end{aligned} \qquad (11) $$

and $\epsilon_z$ the minimum velocity in z-direction. The explanation of $\epsilon_z$ is further detailed in Section 4.3.3.

**Proof.** See Appendix C. ☐

In the following theorems, asymptotic results are given for the average long run travel distance.

**Theorem 2.** The following strong convergences hold for average travel distance,

$$ A_{N_t} \underset{t \to \infty}{\stackrel{\text{a.s.}}{\longrightarrow}} 0, \quad A_t \underset{t \to \infty}{\stackrel{\text{a.s.}}{\longrightarrow}} 0. \qquad (12) $$

**Proof.** See Appendix D. ☐

**Theorem 3.** The following weak convergences hold for average travel distance,

$$ \sqrt{N_t} A_{N_t} \underset{t \to \infty}{\stackrel{D}{\longrightarrow}} \mathcal{N}(0, \sigma_a^2), $$

$$ \sqrt{m_s} A_t \underset{t \to \infty}{\stackrel{D}{\longrightarrow}} \mathcal{N}(0, \sigma_a^2), $$

where

* Under High Temperature

$$ \sigma_a^2 = \frac{1}{2} \mathbb{E}(X_1^2) \theta^2 + \frac{2 - \alpha_{\text{tan}}}{4} \mathbb{E}(X_2^2) \theta^2 + \frac{2m_1^2 + (\alpha_{\text{tan}}^2 - 3\alpha_{\text{tan}} + 4)m_1m_2 + (2 - \alpha_{\text{tan}})m_2^2}{2\alpha_{\text{tan}}} \theta^2, \qquad (13) $$

* Under Low Temperatures

$$ \begin{aligned} \sigma_a^2 = \frac{1}{2} \mathbb{E}(X_1^2) \theta^2 &+ \frac{m_1(m_2 + m_D)}{1 - p(1 - \alpha_{\text{tan}})} \theta^2 + \frac{2 - \alpha_{\text{tan}}}{4} \mathbb{E}(X_2^2 | X_2 < t_0) \theta^2 \\ &+ \frac{1}{2} \mathbb{E}(X_2 | X_2 < t_0) \frac{p(1 - \alpha_{\text{tan}})m_D + m_2}{1 - p(1 - \alpha_{\text{tan}})} \theta^2 + \frac{1}{4} q t_0^2 \theta^2 + \frac{q D_s}{\lambda} + \frac{1}{2} q t_0 \theta^2, \end{aligned} \qquad (14) $$
with

$$m_D = \frac{h}{\theta} \sqrt{\frac{\pi}{2} (1 - erf(\frac{h}{\sqrt{2\theta t_0})))} + (1 - e^{-\frac{h^2}{2\theta^2 t_0^2}}) t_0,$$

$$\mathbb{E}(X_1^2) = -\frac{h^2}{2\theta^2} Ei(-\frac{e_z^2}{2\theta^2}),$$

$$\mathbb{E}(X_2^2) = \frac{p}{e^{-(t_0^2 - m^2)/\theta^2}} \int_0^{t_0} t^2 dF_0(t) + q(t_0^2 + \frac{2t_0}{\lambda} + \frac{2}{\lambda^2}),$$

$$\mathbb{E}(X_2|X_2 < t_0) = \frac{p}{e^{-(t_0^2 - m^2)/\theta^2}} \int_0^{t_0} t dF_0(t),$$

$$\mathbb{E}(X_2^2|X_2 < t_0) = \frac{p}{e^{-(t_0^2 - m^2)/\theta^2}} \int_0^{t_0} t^2 dF_0(t),$$

and $\epsilon_z$ the minimum velocity in z-direction.

**Proof.** See Appendix E. ☐

**Remark.** Using Theorem 3, we can deduce the asymptotic confidence interval for the average travel distance $A_t$ with a confidence level $1 - \eta$, which is

$$[-\sqrt{2\sigma_a erf^{-1}(1-\eta)} \quad \sqrt{2\sigma_a erf^{-1}(1-\eta)}]. \qquad (15)$$

We have $1 - \eta$ percentage of confidence that the particle will be in this interval. In other words, among all particles in the system, there is a proportion $1 - \eta$ of particles which fall within this confidence interval distance.

**Theorem 4.** The diffusion coefficient of channel is given by:

$$D_{\text{channel}} = \frac{\sigma_a^2}{2m_s}. \qquad (16)$$

**Proof.** See Appendix F. ☐

# 4. Numerical application

In this section, we aim to validate the analytical results presented in Section 3. Since the results are categorized into two cases — high temperature and low temperature — we begin by MC simulations for both scenario and comparing their respective asymptotic outcomes. However, as indicated in [33], significant discrepancies arise between the MC simulation and MD results, with deviations reaching approximately 20% at high temperatures and increasing further at low temperatures. At lower temperatures, the diffusion process becomes more pronounced, causing the second-order Fréchet distribution shape parameter to lose validity. Moreover, the transition time $t_0$ is selected empirically rather than derived through a rigorous scientific framework. We believe that the deviation between MC simulation and MD results is partially due to their parameter estimation. Therefore, a more precise estimation of the parameters for the distribution function $F_T$ is important. Here, we present the most effective method along with its verification for this estimation, while alternative approaches are detailed in the Appendix G for completeness. Besides, the proposed method is applied to the MD data to estimate the model parameters more accurately, followed by the calculation of all related analytical theoretical results using the estimated parameters.

## 4.1. Numerical comparison with Monte Carlo simulation

In this section the asymptotic theoretical results obtained in the previous section are compared with MC simulation results of the physical model, in order to verify the accuracy of analytical results. The MC simulation is based on the physical model described in Section 2. The process begins in state 1, where the initial velocity distribution in the longitudinal and transversal directions is assigned. The residence time in state 1 is determined by the particle's velocity in the longitudinal direction, specifically $h/V_z$. Upon entering state 2, the particle's velocity in the x-direction is updated according to the relation in Eq. (1), while the velocity in the z-direction is sampled from a Rayleigh distribution. The residence time in state 2 follows the distribution function given in Eq. (4). This sequence is repeated $N$ times, where $N$ is chosen to be sufficiently large to ensure statistical convergence. The parameters are estimated through regression analysis from [33], with the Fréchet distribution set to $s = 2$ and $m = 0$. The analytical results are presented under two situations in Section 3. Accordingly, 350 K, representing a high-temperature regime where the system exhibits predominantly ballistic behavior, and 50 K, representing a low-temperature regime where both ballistic and diffusive processes coexist. The model parameters are used to generate data by MC simulations and their asymptotic behavior are compared with the analytical results presented in Section 3.
<figure><img src="image_4.png" /><figcaption>Fig. 3. 350 K resident time average $m_s$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. 350 K resident time variance $\sigma^2$.</figcaption></figure>

<table><caption>Table 1
Parameters for 350 K.</caption><thead><tr><th>h (Å)</th><th>z<sub>cut</sub> (Å)</th><th>t<sub>0</sub> (ps)</th><th>α<sub>tan</sub></th><th>β (ps)</th><th>1/λ (ps<sup>-1</sup>)</th><th>p</th></tr></thead><tbody><tr><td>50</td><td>12</td><td>15</td><td>0.045</td><td>1.698</td><td>4.59</td><td>0.985</td></tr></tbody></table>

*High temperature.* For a temperature of 350 K, the parameter settings are listed in Table 1. These parameters are estimated by fitting a regression model on data as in [33]. The theoretical results corresponding to these parameters are presented in Table 2. At a confidence level of 1 − η, the confidence interval for the average travel distance, as described in Theorem 3 is given in Table 3. MC results for the high temperature are presented in Fig. 3, 4, 5, 6, where the convergence of the relevant quantities can be observed. The limiting parameters are summarized also in Table 2 for comparison. Besides, to verify Lemma 2, we compare also 1000 particles’ trajectory $\sqrt{N_t}A_{N_t}$ and the theoretical pdf $N(0, \sigma_a^2)$ in Fig. 7. The analytical asymptotic results have minor relative errors compared to empirical estimates obtained via Monte Carlo simulations. Since variance is the square of the standard deviation, the relative error for the standard deviation is also expected to be small. Based on the relative error results, we reach a similar conclusion that the asymptotic analytical results align very well with the empirical estimates obtained from Monte Carlo simulations.
<figure><img src="image_6.png" /><figcaption>Fig. 5. 350 K travel distance average $A_t$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. 350 K travel distance variance $\sigma_a^2$.</figcaption></figure>

**Table 2**

Analytical results given in Eq. (8), (11), (12), (13) and empirical estimation results considering Monte Carlo simulation from Figs. 3, 4, 5, 6 for high temperature 350 K considering parameters given in Table 1.

<table><thead><tr><th></th><th>m<sub>s</sub> (ps)</th><th>σ<sup>2</sup> (ps<sup>2</sup>)</th><th>ā (Å)</th><th>σ<sub>a</sub><sup>2</sup> (Å<sup>2</sup>)</th></tr></thead><tbody><tr><td>Analytical asymptotic results</td><td>3.3675</td><td>33.0315</td><td>0</td><td>75042.3753</td></tr><tr><td>Empirical estimation results</td><td>3.3677</td><td>33.7960</td><td>0.0009</td><td>72499.3227</td></tr><tr><td>Relative error</td><td>0.006%</td><td>2.3%</td><td>0.09%</td><td>3.5%</td></tr></tbody></table>

*Low temperature.* For a temperature of 50 K, the parameter settings are listed in Table 4 from [33]. The analytical results are exposed in Table 5. At a confidence level of $1 - \eta$, the confidence interval for the average travel distance is given in Table 6. MC results for the high temperature are presented in Figs. 8, 9,10, 11. The figures show the stabilization of different quantity of interest in Table 5. To verify Lemma 2 under low temperature, we compare also 1000 particles’ trajectory $\sqrt{N_t}A_{N_t}$ and the theoretical pdf $N(0, \sigma_a^2)$ in Fig. 12.

Therefore, the model can accurately describe the system dynamics with the appropriate parameters. It proves that the model is capable of capturing essential physical properties.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Histogram of $\sqrt{N_t}A_{N_t}$ compare with theoretical pdf for 350 K.</figcaption></figure>

<caption>**Table 3**</caption>

Analytical asymptotic Confidence Intervals of level 1-η for Aᵢ by Eq. (15) at high temperature 350 K.

<table><thead><tr><th>1-η %</th><th>Confidence interval (Å)</th></tr></thead><tbody><tr><td>80</td><td>[-352.9910, 352.9910]</td></tr><tr><td>85</td><td>[-396.5051, 396.5051]</td></tr><tr><td>90</td><td>[-453.0591, 453.0591]</td></tr><tr><td>95</td><td>[-539.8532, 539.8532]</td></tr><tr><td>98</td><td>[-640.7701, 640.7701]</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 8. 50 K resident time average $m_s$.</figcaption></figure>

<caption>**Table 4**</caption>

Parameters used in [33] for 50 K.

<table><thead><tr><th>d (Å)</th><th>z<sub>cut</sub> (Å)</th><th>α<sub>tan</sub></th><th>β (ps)</th><th>1/λ (ps)</th><th>t<sub>0</sub> (ps)</th><th>p</th><th>D<sub>s</sub> (Å<sup>2</sup>/ps)</th></tr></thead><tbody><tr><td>50</td><td>12</td><td>0.23</td><td>4.157</td><td>140</td><td>60</td><td>0.7</td><td>289</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 9. 50 K resident time variance $\sigma^2$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. 50K travel distance average $A_t$.</figcaption></figure>

<table><caption>**Table 5**
Analytical results given in Eqs. (8), (11), (12), (14) for low temperature 50 K considering parameters given in Table 4.</caption><thead><tr><th></th><th><i>m</i><sub>z</sub> (ps)</th><th>σ<sup>2</sup> (ps<sup>2</sup>)</th><th><i>a</i> (Å)</th><th>σ<sup>2</sup><sub>a</sub> (Å<sup>2</sup>)</th></tr></thead><tbody><tr><td>Analytical asymptotic results</td><td>37.4428</td><td>7861.4063</td><td>0</td><td>38805.3923</td></tr><tr><td>Empirical estimation results</td><td>37.4415</td><td>7971.7521</td><td>0.00006</td><td>38018.5196</td></tr><tr><td>Relative error</td><td>0.003%</td><td>1.4%</td><td>0.006%</td><td>2.1%</td></tr></tbody></table>

## 4.2. Time threshold detection and parameter estimation methods

Although the analytical results correspond to MC simulations, there is still a difference between MC simulation results and MD simulation results according to [33]. The transition point $t_0$ is a key parameter in the analytical results. However, in [33], it is selected empirically rather than determined through a systematic or scientific method. The remaining parameters are then estimated via regression once $t_0$ has been chosen. Therefore, here we present a method to estimate $t_0$ along with other parameters in distribution function $F_T$. There are some alternative methods to this end, presented in Appendix G, and their verification along with comparison is in Appendix H.
<figure><img src="image_12.png" /><figcaption>Fig. 11. 50 K travel distance variance $\sigma_a^2$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Histogram of $\sqrt{N_t A_N}$, compare with theoretical pdf for 50 K.</figcaption></figure>

<caption>Table 6</caption>

Analytical asymptotic Confidence Intervals of level 1-η for A, by Eq. (15) at low temperature 50 K.

<table><thead><tr><th>1 - η (%)</th><th>Confidence interval (Å)</th></tr></thead><tbody><tr><td>80</td><td>[-252.4539, 252.4539]</td></tr><tr><td>85</td><td>[-283.5745, 283.5745]</td></tr><tr><td>90</td><td>[-324.0211, 324.0211]</td></tr><tr><td>95</td><td>[-386.0950, 386.0950]</td></tr><tr><td>98</td><td>[-458.2692, 458.2692]</td></tr></tbody></table>

*Empirical cumulative distribution function.* Considering the whole data, the analytical cdf corresponding to the whole model $F_T$ by Eq. (4) is compared with the empirical cumulative distribution $\hat{F}_T$. Since the cdf depends on the parameters and the time threshold
<table><caption>**Table 7**
Imposed parameter values for Monte Carlo simulation and analytical theoretical results.</caption><thead><tr><th></th><th>s</th><th>β (ps)</th><th>m (ps)</th><th>λ (ps)</th><th>p</th><th>t<sub>0</sub> (ps)</th></tr></thead><tbody><tr><td>Parameters for different mean</td><td>1.5</td><td>2.0</td><td>1.0</td><td>0.05</td><td>0.6</td><td>20</td></tr><tr><td>Parameters for similar mean</td><td>1.5</td><td>2.0</td><td>1.0</td><td>0.15</td><td>0.6</td><td>20</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Error with different time threshold with different mean by MLE (left) by minimization of cdf (right).</figcaption></figure>

<table><caption>**Table 8**
Error with different time threshold with different mean by MLE and by minimization of cdf.</caption><thead><tr><th></th><th>t̂<sub>0</sub></th><th>ŝ</th><th>β̂ (ps)</th><th>m̂ (ps)</th><th>λ̂ (1/ps)</th><th>β̂</th></tr></thead><tbody><tr><td>MLE</td><td>19.6970</td><td>1.4866</td><td>1.9981</td><td>1.0041</td><td>0.0492</td><td>0.5985</td></tr><tr><td>Minimization of cdf</td><td>20.0404</td><td>1.5036</td><td>2.0002</td><td>0.9970</td><td>0.0502</td><td>0.5998</td></tr></tbody></table>

$t_0$, all the unknown quantities can be estimated by minimizing over the parameters and $t_0$ the following distance

$$E = \sum_{i=1}^{N} (\hat{F}_{T}(t_i) - F_{T}(t_i))^2, \quad \text{where} \quad \hat{F}_{T}(t) = \frac{1}{N} \sum_{i=1}^{N} \mathbb{1}(t_i \le t). \qquad (17)$$

Since the probability $p$ is the weight permitting the normalization, once the transition time $t_0$ is determined, it can be estimated as follows:

$$\hat{p} = \frac{1}{N} \sum_{i=1}^{N} \mathbb{1}(t_i \le t_0).$$

The remaining parameters, namely $s$, $\beta$, and $m$ for the Fréchet distribution, and $\lambda$ for the Exponential distribution, are estimated using either the Maximum Likelihood Estimation (MLE) method or by minimizing the mean square error of the cumulative distribution function estimator.

In order to verify the estimation methods, we apply these latter to MC simulation data and compare the estimated parameters with values. We generate $N = 100,000$ collisions by MC Simulations in the same way presented in Section 4. We propose two sets of imposed values for parameters in Table 7. The chosen values are within the same order of magnitude as the estimated parameters listed in Table 8. To better observe the transition dynamics and appropriately weight the Fréchet and Exponential components, slight adjustments were made to the parameter values. Besides, it is generally known in statistics that the closer the means before and after the threshold, the more difficult it becomes to estimate the change time. Therefore, we apply the same method to different cases where the two distributions have close mean values, in order to justify different situations. These generated data are used for time threshold estimation and parameter estimation.

In this method, we estimate the 4 parameters of $F_T$ with MLE or with minimization of cdf error given in Eq. (17).

The distance $E$ is depicted and numerically we get the minimum at 19.6970 by MLE and at 20.0404 of cdf error given in Eq. (17). There is a short gap between results and they are close to our real value 20 usd for simulation. These methods have been applied a large number of times and the results are as promising as what is exposed in Fig. 13. (See Table 8.)

It can be pointed out that in the previous case the mean for Fréchet distribution and the mean for Exponential distribution are very different. For close means, the results are depicted in Fig. 14, where the minimum of cdf distance is obtained at 20 ps when the parameters estimation is carried out by MLE (left) which has the minimum at 20.026 ps. The parameters estimation is carried out with minimization of cdf in Fig. 14 (right) which has the minimum at 10 ps. Even if distance minimization does not permit the detection time, the detection is possible. These methods have been applied a large number of times and the results are as promising as what is exposed in the Figures when the parameters estimation is carried out by MLE.
<figure><img src="image_15.png" /><figcaption>Fig. 14. Error with different time threshold with similar mean by MLE (left) by minimization of cdf (right).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Estimated cdf and empirical cdf comparison for 350 K (left) and for 50 K (right).</figcaption></figure>

This method provides the best estimation with precision for time threshold, and MLE performs well in all situation whether the mean of two distributions are close. Other methods in Appendix G give vaguely an interval, but lack sufficient precision.

### 4.3. Molecular dynamics data

We now apply the proposed estimation method to MD simulation results to evaluate whether it improves the accuracy of the analytical predictions. We generated a series of data through MD simulations for temperatures of 350 K, 200 K for high temperature, and 50 K for low temperature. The 200 K condition is categorized as high temperature since few particles exhibit diffusion processes, making it excessive to treat it as low temperature. As described in the Introduction, these simulations focus on the He/graphite system under thermal equilibrium, with gas/gas interactions neglected (infinite Knudsen number). There are 40 Helium particles in the channel which generate 200000 collisions. Each collision has a resident time, and the data is formed in group data. From 0 to 250 ps, the resident time is divided into small time interval of 0.125 ps. MD gives the number of particles with corresponding time interval.

#### 4.3.1. Estimation for molecular dynamics data

In this subsection, we present the estimation process for high temperature 350 K and low temperature 50 K. The intermediate temperature of 200 K is treated as high temperature, as the diffusion behavior of several particles aligns more closely with high-temperature dynamics. Since MD data are interval-censored, the likelihood is written as follows [47–49]:

$$ \mathcal{L} = \prod_{i=1}^{N-1} (\hat{F}(x_i) - \hat{F}(x_{i-1}))^{N_i} (\bar{F}(x_N))^{N_N}, $$

where $N_i$ is the number of particles in time interval $[x_{i-1}, x_i]$.

By maximizing the likelihood $\mathcal{L}$, the model parameters with generalized Fréchet distribution are estimated.

The estimated cdf and empirical cdf for high temperature 350 K and low temperature 50 K is compared in Fig. 15.
<table><caption>**Table 9**
Estimation results for MD data considering interval-censored data.</caption><thead><tr><th>Temperature (K)</th><th>Parameter</th><th>Generalized Fréchet</th><th>Fréchet order 2</th></tr></thead><tbody><tr><td rowspan="6">350</td><td>ŷ<sub>0</sub></td><td>4.0804</td><td>5.0004</td></tr><tr><td>ŝ</td><td>2.0791</td><td>2</td></tr><tr><td>β̂ (ps)</td><td>1.5661</td><td>1.6984</td></tr><tr><td>m̂ (ps)</td><td>0.1943</td><td>0</td></tr><tr><td>λ̂ (1/ps)</td><td>0.2595</td><td>0.2486</td></tr><tr><td>p̂</td><td>0.8044</td><td>0.8308</td></tr><tr><td rowspan="2">200</td><td>α<sub>tan</sub></td><td>0.045</td><td>0.045</td></tr><tr><td>D<sub>s</sub> (Å<sup>2</sup>/ps)</td><td>-</td><td>-</td></tr><tr><td rowspan="6">50</td><td>ŷ<sub>0</sub></td><td>6.1709</td><td>6.3098</td></tr><tr><td>ŝ</td><td>1.6883</td><td>2</td></tr><tr><td>β̂ (ps)</td><td>1.8388</td><td>2.1743</td></tr><tr><td>m̂ (ps)</td><td>0.3158</td><td>0</td></tr><tr><td>λ̂ (1/ps)</td><td>0.1536</td><td>0.1531</td></tr><tr><td>p̂</td><td>0.7800</td><td>0.8520</td></tr><tr><td rowspan="2">200</td><td>α<sub>tan</sub></td><td>0.066</td><td>0.066</td></tr><tr><td>D<sub>s</sub> (Å<sup>2</sup>/ps)</td><td>-</td><td>-</td></tr><tr><td rowspan="6">50</td><td>ŷ<sub>0</sub></td><td>32.1000</td><td>22.0381</td></tr><tr><td>ŝ</td><td>1.1376</td><td>2</td></tr><tr><td>β̂ (ps)</td><td>2.8057</td><td>3.7966</td></tr><tr><td>m̂ (ps)</td><td>1.2721</td><td>0</td></tr><tr><td>λ̂ (1/ps)</td><td>0.0085</td><td>0.0078</td></tr><tr><td>p̂</td><td>0.6347</td><td>0.5986</td></tr><tr><td rowspan="2">200</td><td>α<sub>tan</sub></td><td>0.184</td><td>0.164</td></tr><tr><td>D<sub>s</sub> (Å<sup>2</sup>/ps)</td><td>289</td><td>-</td></tr></tbody></table>

<table><caption>**Table 10**
Analytical results with estimated parameters by MD data compared with MD data.</caption><thead><tr><th>T (K)</th><th></th><th>Estimation with Fréchet order 2</th><th>Estimation with Generalized Fréchet</th><th>MD</th></tr></thead><tbody><tr><td rowspan="2">350</td><td>E(X<sub>2</sub>)</td><td>3.0765</td><td>3.2719</td><td>3.1400</td></tr><tr><td>Var(X<sub>2</sub>)</td><td>19.2778</td><td>10.1297</td><td>10.0533</td></tr><tr><td rowspan="2">200</td><td>E(X<sub>2</sub>)</td><td>5.4071</td><td>4.8637</td><td>4.8838</td></tr><tr><td>Var(X<sub>2</sub>)</td><td>31.8567</td><td>27.6293</td><td>28.6653</td></tr><tr><td rowspan="2">50</td><td>E(X<sub>2</sub>)</td><td>67.8887</td><td>59.1118</td><td>63.5414</td></tr><tr><td>Var(X<sub>2</sub>)</td><td>12336.2622</td><td>9804.1042</td><td>8596.1311</td></tr></tbody></table>

We have also conducted similar estimations using a Fréchet distribution of order 2, with parameters $s = 2$ and $m = 0$, to enable a comparison with the generalized Fréchet distribution and their respective outcomes. The results of these estimations are presented in Table 9. Additionally, the corresponding system parameter, TMAC ($α_t$), and the surface diffusion coefficient ($D_s$) are included in the same table. It should be noticed that the TMAC values used here correspond to the regime before $t_0$, as defined in Eq. (1). Consequently, these values vary depending on changes in $t_0$, as shown in the table. Under high temperature, it is assumed that there is no surface diffusion, therefore no $D_s$ for high temperature.

### 4.3.2. Asymptotic results based in estimated parameters

In this section, we employ empirical cumulative distribution function method for parameter estimation. Therefore, we use the estimation results in Table 9 to calculate the asymptotic mean and variance of resident time and compare them with MD data in Table 10.

We give the asymptotic results in Table 11, the related confident interval is given in Table 12. Here, $ε_z$ is chosen by a minimum of MC simulation of Rayleigh distribution.

Using Eq. (16), the diffusion coefficient of the channel can be calculated analytically, and the results are compared with MD simulations in Table 13. The analytical error remains within an acceptable range. The generalized Fréchet distribution demonstrates better performance at high temperatures. Although its accuracy slightly decreases at low temperatures, it remains within an acceptable range. This behavior is likely due to the influence of the time threshold parameter, $t_0$, which affects the TMAC estimation. In [33], the authors adjust the transition time $t_0$ and propose several diffusion coefficient results. Some of them have the same order of MD results. However, MC simulation also take time to get a result, while analytical expressions allow for efficient computation of results based on system parameters and enable a systematic analysis of the influence of each parameter on the asymptotic outcomes.
<table><caption>**Table 11**

Analytical results with estimated parameters in Table 9.</caption><thead><tr><th>Condition</th><th></th><th>m<sub>s</sub> (ps)</th><th>σ<sup>2</sup> (ps<sup>2</sup>)</th><th>ā (Å)</th><th>σ<sub>a</sub><sup>2</sup> (Å<sup>2</sup>)</th></tr></thead><tbody><tr><td rowspan="2">350 K,<br>ε<sub>z</sub> = 0.004 (Å/ps)</td><td>Fréchet order 2</td><td>3.5462</td><td>33.7539</td><td>0</td><td>84226.7268</td></tr><tr><td>Generalized Fréchet</td><td>3.5528</td><td>33.7711</td><td>0</td><td>83244.2273</td></tr><tr><td rowspan="2">200 K,<br>ε<sub>z</sub> = 0.004 (Å/ps)</td><td>Fréchet order 2</td><td>4.9588</td><td>61.6062</td><td>0</td><td>63512.7817</td></tr><tr><td>Generalized Fréchet</td><td>4.5434</td><td>63.1235</td><td>0</td><td>70402.3514</td></tr><tr><td rowspan="2">50 K,<br>ε<sub>z</sub> = 0.0004 (Å/ps)</td><td>Fréchet order 2</td><td>34.6099</td><td>5745.7173</td><td>0</td><td>31953.9587</td></tr><tr><td>Generalized Fréchet</td><td>38.9983</td><td>7246.1048</td><td>0</td><td>37932.5619</td></tr></tbody></table>

<table><caption>**Table 12**

Analytical asymptotic Confidence Intervals of level 1 − η for A<sub>γ</sub>.</caption><thead><tr><th>1 − η (%)</th><th>Model</th><th>Temperature</th><th>Confidence Interval (Å)</th></tr></thead><tbody><tr><td rowspan="6">80</td><td rowspan="3">Fréchet order 2</td><td>350 K</td><td>[−372.9299, 372.9299]</td></tr><tr><td>200 K</td><td>[−322.9733, 322.9733]</td></tr><tr><td>50 K</td><td>[−229.0859, 229.0859]</td></tr><tr><td rowspan="3">Generalized Fréchet</td><td>350 K</td><td>[−369.7542, 369.7542]</td></tr><tr><td>200 K</td><td>[−340.0397, 340.0397]</td></tr><tr><td>50 K</td><td>[−249.5986, 249.5986]</td></tr><tr><td rowspan="6">85</td><td rowspan="3">Fréchet order 2</td><td>350 K</td><td>[−417.7785, 417.7785]</td></tr><tr><td>200 K</td><td>[−362.7870, 362.7870]</td></tr><tr><td>50 K</td><td>[−257.3259, 257.3259]</td></tr><tr><td rowspan="3">Generalized Fréchet</td><td>350 K</td><td>[−415.3347, 415.3347]</td></tr><tr><td>200 K</td><td>[−381.9572, 381.9572]</td></tr><tr><td>50 K</td><td>[−280.3672, 280.3672]</td></tr><tr><td rowspan="6">90</td><td rowspan="3">Fréchet order 2</td><td>350 K</td><td>[−477.3368, 477.3368]</td></tr><tr><td>200 K</td><td>[−414.5318, 414.5318]</td></tr><tr><td>50 K</td><td>[−294.0286, 294.0286]</td></tr><tr><td rowspan="3">Generalized Fréchet</td><td>350 K</td><td>[−474.5744, 474.5744]</td></tr><tr><td>200 K</td><td>[−436.4363, 436.4363]</td></tr><tr><td>50 K</td><td>[−320.3564, 320.3564]</td></tr><tr><td rowspan="6">95</td><td rowspan="3">Fréchet order 2</td><td>350 K</td><td>[−586.8176, 586.8176]</td></tr><tr><td>200 K</td><td>[−493.9451, 493.9451]</td></tr><tr><td>50 K</td><td>[−350.3567, 350.3567]</td></tr><tr><td rowspan="3">Generalized Fréchet</td><td>350 K</td><td>[−565.4903, 565.4903]</td></tr><tr><td>200 K</td><td>[−520.1459, 520.1459]</td></tr><tr><td>50 K</td><td>[−381.7281, 381.7281]</td></tr><tr><td rowspan="6">98</td><td rowspan="3">Fréchet order 2</td><td>350 K</td><td>[−675.1490, 675.1490]</td></tr><tr><td>200 K</td><td>[−586.2802, 586.2802]</td></tr><tr><td>50 K</td><td>[−415.8503, 415.8503]</td></tr><tr><td rowspan="3">Generalized Fréchet</td><td>350 K</td><td>[−671.1996, 671.1996]</td></tr><tr><td>200 K</td><td>[−617.2602, 617.2602]</td></tr><tr><td>50 K</td><td>[−453.0861, 453.0861]</td></tr></tbody></table>

<table><caption>**Table 13**

Analytical diffusion coefficient of channel D<sub>channel</sub> compared with MD simulation results.</caption><thead><tr><th></th><th></th><th>MD D<sub>channel</sub></th><th>Analytical D<sub>channel</sub></th><th>Relative error</th></tr></thead><tbody><tr><td rowspan="2">350 K</td><td>Fréchet order 2</td><td rowspan="2">12640</td><td>11443.6604</td><td>9.5%</td></tr><tr><td>Generalized Fréchet</td><td>12352.8054</td><td>2.3%</td></tr><tr><td rowspan="2">200 K</td><td>Fréchet order 2</td><td rowspan="2">6431</td><td>5929.7392</td><td>7.8%</td></tr><tr><td>Generalized Fréchet</td><td>6404.0475</td><td>0.4%</td></tr><tr><td rowspan="2">50 K</td><td>Fréchet order 2</td><td rowspan="2">480</td><td>486.3361</td><td>1.3%</td></tr><tr><td>Generalized Fréchet</td><td>497.2074</td><td>3.6%</td></tr></tbody></table>

### 4.3.3. Sensitivity analyses of ε<sub>z</sub>

In the above calculation, we use the minimum v<sub>z</sub> simulation result to be ε<sub>z</sub>. To address the influence of ε<sub>z</sub> on results, we perform a sensitivity analysis with 75 K. Its influence on σ<sup>2</sup> and on σ<sub>a</sub><sup>2</sup> is in Fig. 16.

As observed by Arya et al. [35] a transversal velocity cutoff must be introduced to compute the mean square time spent in state 1. Arya et al. found also a logarithmic cutoff dependence of mean square time between the two wall collisions and of the self-diffusion. In the case of a cylindrical and a rectangular channel of finite length, Smoluchowski [34] computed the gas flux taking into account of the TMAC. In the case of the rectangular geometry, Smoluchowski encountered also a logarithmic divergence of the diffusion coefficient as the side ratio tends to infinity. The figures show that as ε<sub>z</sub> increases, both σ<sup>2</sup> and σ<sub>a</sub><sup>2</sup> decrease. However, the results
<figure><img src="image_17.png" /><figcaption>Fig. 16. Sensitivity of $ε_z$ on resident time variance $σ^2$ (left) and on travel distance variance $σ_u^2$ (right).</figcaption></figure>

remain in the same order, so that the influence of $ε_z$ appears to be limited as observed by Arya et al. So, the slit pore is not the best geometry to investigate the gas diffusion process especially as the gas/gas interaction is neglected, but it allows to access to very simple geometry.

## 5. Conclusion

In this paper, we proposed a semi-Markov model to describe the movement of helium particles in a graphite channel, considering both high and low temperature regimes. The numerical results obtained from the theoretical model closely align with those produced by the MC simulation for physical model. These analytical formulas not only facilitate rapid computations based on system parameters but also enable a comprehensive understanding of how each parameter affects the asymptotic results, offering valuable insights into the dynamics of the system.

Given the discrepancy between the MC and MD simulations, and the critical role of the transition time $t_0$ in the analytical theoretical results, we propose several methods to enhance parameter estimation in order to minimize this discrepancy. These methods are rigorously validated against MC simulation results. Besides, the most effective method is employed to estimate the parameters for the MD simulation. Then the corresponding analytical theoretical results are computed using the estimated parameters.

Regarding the parameter $ε_z$, which represents the minimum velocity in the z-direction, which permits a particle to enter the interaction zone, we study the sensitivity of $ε_z$ and find that while it influences the results, the asymptotic outcomes remain of the same order, indicating a limited impact.

For future work, we aim to explore two main directions: further refinement of the theoretical model, including convergence rates and deviation results, as well as more robust statistical estimation methods to identify parameters that better correspond to real-world conditions. On the physical modeling side, extending the analysis to calculate the diffusion coefficient is a promising avenue. Additionally, considering the transition between ballistic and diffusive regimes, as well as the potential inclusion of sub-diffusion processes, could be fruitful for future studies. We also plan to extend the system's complexity by considering a cylindrical geometry of other material and accounting for gas-gas collisions, which would enhance the model's realism in practical applications.

### CRediT authorship contribution statement

**Jingqi Zhang:** Writing – review & editing, Writing – original draft, Validation, Supervision, Methodology, Investigation, Formal analysis, Data curation. **Mitra Fouladirad:** Writing – review & editing, Validation, Supervision, Project administration. **Nikolaos Limnios:** Writing – review & editing, Validation, Supervision, Methodology, Investigation, Formal analysis, Conceptualization. **Pierre Magnico:** Writing – review & editing, Validation, Supervision, Software, Formal analysis, Data curation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgment

The authors are grateful to the reviewers for their valuable comments, which have helped improve the presentation of this work. This work is partially supported by the China Scholarship Council. The “Centre de Calcul Intensif d’Aix-Marseille” is acknowledged for granting access to its high performance computing resources.
# Appendix A. Proof of Proposition 1

With function (1), we can calculate the average and variance of $V_x$ directly:

$$ 
\begin{aligned}
\mathbb{E}(V_x^{(n)}) &= (1 - \alpha_{\text{tan}})\mathbb{E}(V_x^{(n-1)}) + (\alpha_{\text{tan}}(2 - \alpha_{\text{tan}}))^{\frac{1}{2}}\mathbb{E}(V_x^M) \\
&= (1 - \alpha_{\text{tan}})^2\mathbb{E}(V_x^{(n-2)}) + (\alpha_{\text{tan}}(2 - \alpha_{\text{tan}}))^{\frac{2}{2}}\mathbb{E}(V_x^M) \\
&= \dots \\
&= (1 - \alpha_{\text{tan}})^n\mathbb{E}(V_x^{(0)}) + (\alpha_{\text{tan}}(2 - \alpha_{\text{tan}}))^{\frac{n}{2}}\mathbb{E}(V_x^M) = 0,
\end{aligned}
$$

$$ 
\begin{aligned}
\text{Var}(V_x^{(n)}) &= (1 - \alpha_{\text{tan}})^2 \text{Var}(V_x^{(n-1)}) + \alpha_{\text{tan}}(2 - \alpha_{\text{tan}}) \text{Var}(V_x^M) \\
&= (1 - \alpha_{\text{tan}})^4 \text{Var}(V_x^{(0)}) + ((1 - \alpha_{\text{tan}})^2 + 1) \alpha_{\text{tan}}(2 - \alpha_{\text{tan}}) \text{Var}(V_x^M) \\
&= \dots \\
&= (1 - \alpha_{\text{tan}})^{2n} \text{Var}(V_x^{(0)}) + ((1 - \alpha_{\text{tan}})^{2n-2} + \dots + 1) (\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})) \text{Var}(V_x^M) \\
&= \theta^2.
\end{aligned}
$$

# Appendix B. Proof of Theorem 1

The following Lemma is used in this proof.

**Lemma 2.** If $X_n \sim \mathcal{N}(\mu_n, \sigma_n^2)$, $\mu_n \in \mathbb{R}$, and $\sigma_n^2 > 0$ with $\mu_n \to \mu$, $\sigma_n^2 \to \sigma^2$, then we have

$$ X_n \xrightarrow{D} X, \quad \text{where} \quad X \sim \mathcal{N}(\mu, \sigma^2). $$

**Stationary law solution:**

To find the explicit expression of $\rho$, we take $B = (-\infty, x)$ of Eq. (9).

(a) Under high temperature,

$$ \int_{-\infty}^{x} \rho(u) du = \int_{-\infty}^{+\infty} \rho(u) du \int_{-\infty}^{x} f_N \left( \frac{v - (1 - \alpha_{\text{tan}})u}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} \right) \frac{dv}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} $$

Deriving from both sides, we obtain

$$ \rho(x) = \frac{1}{\theta \sqrt{2\pi} \alpha_{\text{tan}}(2 - \alpha_{\text{tan}})} \int_{-\infty}^{+\infty} \rho(u) e^{-\frac{(x - (1 - \alpha_{\text{tan}})u)^2}{\theta \sqrt{2\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}} du. \quad (B.1) $$

A solution of Eq. (B.1) is

$$ \rho(x) = \frac{1}{\theta \sqrt{2\pi}} e^{-\frac{1}{2} \frac{x^2}{\theta^2}} \sim \mathcal{N}(0, \theta^2). $$

(b) Under low temperature

$$ \begin{cases} \int_{-\infty}^{x} \rho(u) du = \int_{-\infty}^{+\infty} \rho(u) du \int_{-\infty}^{x} f_N \left( \frac{v - (1 - \alpha_{\text{tan}})u}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} \right) \frac{dv}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} , & t < t_0, \\ \int_{-\infty}^{x} \rho(u) du = \int_{-\infty}^{+\infty} \rho(u) du \int_{-\infty}^{x} f_N(v) dv, & t \ge t_0. \end{cases} $$

Deriving from both sides, we have:

$$ \begin{cases} \rho(x) = \frac{1}{\theta \sqrt{2\pi} \alpha_{\text{tan}}(2 - \alpha_{\text{tan}})} \int_{-\infty}^{+\infty} \rho(u) e^{-\frac{(x - (1 - \alpha_{\text{tan}})u)^2}{\theta \sqrt{2\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}} du, & t < t_0, \\ \rho(x) = \frac{1}{\theta \sqrt{2\pi}} e^{-\frac{1}{2} \frac{x^2}{\theta^2}} \int_{-\infty}^{+\infty} \rho(u) du, & t \ge t_0. \end{cases} \quad (B.2) $$

Then we obtain the result:

$$ \rho(x) = \frac{1}{\theta \sqrt{2\pi}} e^{-\frac{1}{2} \frac{x^2}{\theta^2}} \sim \mathcal{N}(0, \theta^2). \quad (B.3) $$

Therefore, Eq. (B.3) is valid for all temperature.

**Unique solution:**

Since with Eq. (5), the EMC is ergodic. And with Eq. (8), $\mathbb{E}(X_n) < \infty$. Besides, Proposition 1 gives that $\mathbb{E}(V_x) < \infty$. Therefore, the extended SMP is also ergodic, which means solution $\rho$ is the unique solution.
## Appendix C. Proof of Lemma 1

With the formula of total variance [50], we have:

$$ 
\begin{aligned}
\sigma^2 ={}& \mathrm{Var}(X_n) = \mathrm{Var}(X_n|J_n=1)\mathbb{P}(J_n=1) + \mathrm{Var}(X_n|J_n=2)\mathbb{P}(J_n=2) \\
& + (\mathbb{E}(X_n|J_n=1))^2(1-\mathbb{P}(J_n=1))\mathbb{P}(J_n=1) \\
& + (\mathbb{E}(X_n|J_n=2))^2(1-\mathbb{P}(J_n=1))\mathbb{P}(J_n=2) \\
& - 2\mathbb{E}(X_n|J_n=1)\mathbb{P}(J_n=1)\mathbb{E}(X_n|J_n=2)\mathbb{P}(J_n=2).
\end{aligned}
$$

Because of the independence of $X_n$, we have:

$$ 
\begin{aligned}
\sigma^2 &= \frac{1}{2}(\mathbb{E}(X_1^2) - \mathbb{E}^2(X_1)) + \frac{1}{2}(\mathbb{E}(X_2^2) - \mathbb{E}^2(X_2)) + \frac{1}{4}\mathbb{E}^2(X_1) \\
&\quad + \frac{1}{4}\mathbb{E}^2(X_2) - \frac{1}{2}\mathbb{E}(X_1)\mathbb{E}(X_2) \\
&= \frac{1}{2}\mathbb{E}^2(X_1) + \frac{1}{2}\mathbb{E}^2(X_2) - \frac{1}{4}\mathbb{E}^2(X_1) - \frac{1}{4}\mathbb{E}^2(X_2) - \frac{1}{2}\mathbb{E}(X_1)\mathbb{E}(X_2).
\end{aligned}
$$

We calculate separately each term, for state 1,

$$ 
\begin{aligned}
\mathbb{E}(X_1) &= m_1 = \frac{h}{\theta}\sqrt{\frac{\pi}{2}}, \\
\mathbb{E}(X_1^2) &= \int_0^\infty t^2 dG(t) = \int_0^\infty t^2 d\left(e^{-\frac{h^2}{2\theta^2 t^2}}\right) \\
&= \int_0^\infty t^2 \frac{h^2}{\theta^2 t^3} e^{-\frac{h^2}{2\theta^2 t^2}} dt = \int_0^\infty \frac{h^2}{\theta^2 t} e^{-\frac{h^2}{2\theta^2 t^2}} dt.
\end{aligned}
$$

In the physical model, we assume that every particle is surly to hit the wall, so that it exists a minimum $z$-direction velocity $c_z$ as the transversal cutoff velocity. As the resident time at state 1 is $h/V_z$, so that the resident time at 1 is smaller than $h/c_z$. Therefore:

$$ 
\begin{aligned}
\mathbb{E}(X_1^2) &= \int_0^{\frac{h}{c_z}} t^2 dG(t) = \int_0^{\frac{h}{c_z}} t^2 d\left(e^{-\frac{h^2}{2\theta^2 t^2}}\right) = \int_0^{\frac{c_z}{h}} t^2 \frac{h^2}{\theta^2 t^3} e^{-\frac{h^2}{2\theta^2 t^2}} dt \\
&= \frac{h^2}{2\theta^2} E_1\left(\frac{c_z^2}{2\theta^2}\right) = -\frac{h^2}{2\theta^2} E_1\left(-\frac{c_z^2}{2\theta^2}\right).
\end{aligned}
$$

For state 2, we have

$$ 
\begin{aligned}
\mathbb{E}(X_2) &= m_2 = \frac{p}{e^{-\left(\frac{t_0-m}{\bar{p}}\right)^{-s}} - s} \int_0^{t_0} t dF_0(t) + q(t_0 + \frac{1}{\lambda}). \\
\mathbb{E}(X_2^2) &= \frac{p}{e^{-\left(\frac{t_0-m}{\bar{p}}\right)^{-s}} - s} \int_0^{t_0} t^2 dF_0(t) + \int_{t_0}^{\infty} t^2 dF_1(t) \\
&= \frac{p}{e^{-\left(\frac{t_0-m}{\bar{p}}\right)^{-s}} - s} \int_0^{t_0} t^2 dF_0(t) + q(t_0^2 + \frac{2t_0}{\lambda} + \frac{2}{\lambda^2}).
\end{aligned}
$$

With all these calculations, we obtain

$$ 
\begin{aligned}
\sigma^2 &= \frac{1}{2}(\mathbb{E}(X_1^2) + \mathbb{E}(X_2^2)) - \frac{1}{4}(m_1^2 + m_2^2 + 2m_1m_2) \\
&= \frac{1}{2}(\mathbb{E}(X_1^2) + \mathbb{E}(X_2^2)) - \frac{1}{4}(m_1 + m_2)^2.
\end{aligned}
$$

## Appendix D. Proof of Theorem 2

As we have $N_t \to \infty$, when $t \to \infty$, the SLLN can be stated as

$$ A_{N_t} = \frac{1}{N_t} \sum_{n=1}^{N_t} a(J_n, X_n, V_x^{(n-1)}, V_x^{(n)}) \xrightarrow{a.s.} \bar{a}. $$
We have:

$$
\begin{align*}
\bar{a} &:= \sum_{i \in E} v_i \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{\pi}(dt, du, dv) a_i(t, u, v) \\
&= \frac{1}{2} \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} G(dt) \rho(u) P_1(u, dv) a_1(t, u, v) du \\
&\quad + \frac{1}{2} \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} F_T(dt) \rho(u) P_2(u, dv) a_2(t, u, v) du \\
&= \frac{1}{2} \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} G(dt) \mathbb{1}_u(dv) \rho(u) ut du \\
&\quad + \frac{1}{2} \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} F_T(dt) \rho(u) f_N \left( \frac{v - (1 - \alpha_{\tan})u}{\sqrt{\alpha_{\tan}(2 - \alpha_{\tan})}} \right) \left( \frac{u + v}{2} t \right) du.
\end{align*}
$$

Because $\rho(u)$ and $P_t(u, dv)$ are even functions, we have $\int_{\mathbb{R}} u\rho(u)du = 0$, $\int_{\mathbb{R}} vP_2(u, dv) = 0$. As rather under high temperature or low temperature, $a$ is a function of $u$ and $v$, besides, the mean of $Y$ is 0, therefore, $\bar{a} = 0$, and

$$A_{N_t} \xrightarrow{a.s.} \bar{a} = 0.$$

For $A_t$, we can change it into the following form,

$$A_t = \frac{1}{t} \sum_{n=1}^{N_t} a(J_n, X_n, V_x^{(n-1)}, V_x^{(n)}) + \frac{a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t)}, V_x^{(N_t+1)})}{t}.$$

By Eq. (10), we have,

$$\frac{N_t}{t} \frac{1}{N_t} \sum_{n=1}^{N_t} a(J_n, X_n, V_x^{(n-1)}, V_x^{(n)}) \xrightarrow{a.s.} \frac{\bar{a}}{m_s} = 0.$$

Let $(t_n)_{n \ge 1}$ be an increasing sequence of positive real numbers, $t_n \uparrow \infty$. By Chebychev's inequality, we have:

$$P\left(\frac{a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)})}{t_k} - \bar{a} \ge \epsilon\right) \le \frac{E\left(a^2(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)})\right)}{t_k^2 \epsilon^2},$$

since $\bar{a} = 0$ aforementioned.

Because,

$$E_i(a^2(J_i, X_i, V_x^{(0)}, V_x^{(1)})) < \infty, \quad i = 1, 2.$$

we have,

$$P\left(\frac{a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)})}{t} \ge \epsilon\right) \le \frac{E\left(a^2(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)})\right)}{t^2 \epsilon^2} \xrightarrow{t \to \infty} 0. \quad (D.1)$$

Then we have,

$$\sum_{k=1}^{\infty} P\left(\frac{a(J_{N_t k+1}, t - S_{N_t k}, V_x^{(N_t k-1)}, V_x^{(N_t k)})}{t_k} \ge \epsilon\right) \le \sum_{k=1}^{\infty} \frac{E\left(a^2(J_{N_t k+1}, t - S_{N_t k}, V_x^{(N_t k-1)}, V_x^{(N_t k)})\right)}{t_k^2 \epsilon^2} \\ \sum_{k=1}^{\infty} P\left(\frac{a(J_{N_t k+1}, t - S_{N_t k}, V_x^{(N_t k-1)}, V_x^{(N_t k)})}{t_k} \ge \epsilon\right) \le \frac{C}{\epsilon^2} \sum_{k=1}^{\infty} \frac{1}{t_k^2} < \infty,$$

$$\text{where } C = \sup_k E(a^2(J_{N_t k+1}, t - S_{N_t k}, V_x^{(N_t k-1)}, V_x^{(N_t k)})) < \infty.$$

Hence, see e.g. in [51], we have,

$$\frac{a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t)}, V_x^{(N_t+1)})}{t} \xrightarrow{a.s.} 0.$$

Therefore, by Slutsky's theorem, we have,

$$A_t \xrightarrow{a.s.} 0.$$

*Remark.* From Chebychev's inequality (D.1), we have a convergence in probability of $a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)})/t$ to 0, which implies that:

$$A_t \xrightarrow{P} 0.$$
## Appendix E. Proof of Theorem 3

By the CLT [36], we have:

$$ \sqrt{N_t}(A_{N_t} - \bar{a}) \xrightarrow{t \to \infty} \mathcal{N}(0, \sigma_a^2). $$

Because $\bar{a} = 0$, so that,

$$ \sqrt{N_t} A_{N_t} \xrightarrow{t \to \infty} \mathcal{N}(0, \sigma_a^2). $$

As for,

$$ \sqrt{m_s t} A_t = \sqrt{m_s t} \frac{1}{t} \sum_{n=1}^{N_t} a(J_n, N_t, V_x^{(n-1)}, V_x^{(n)}) + \sqrt{m_s t} \frac{a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)})}{t}, $$

we have,

$$ \sqrt{m_s t} \frac{1}{t} \sum_{n=1}^{N_t} a(J_n, N_t, V_x^{(n-1)}, V_x^{(n)}) = \sqrt{m_s t} \frac{N_t}{t} A_{N_t} = \sqrt{m_s} \frac{N_t}{t} \sqrt{N_t} A_{N_t}. $$

As stated in Renewal theory, we have,

$$ \sqrt{m_s} \frac{N_t}{t} \xrightarrow{t \to \infty} 1. $$

Therefore

$$ \sqrt{m_s t} \frac{1}{t} \sum_{n=1}^{N_t} a(J_n, N_t, V_x^{(n-1)}, V_x^{(n)}) \xrightarrow{t \to \infty} \mathcal{N}(0, \sigma_a^2). $$

Besides, by Chebychev's inequality,

$$ \mathbb{P}\left(\left|\sqrt{m_s t} \frac{a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)})}{t}\right| \ge \epsilon\right) \le \frac{m_s \mathbb{P}\mathit{var}(a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)}))}{t \epsilon^2}. $$

As $\mathbb{P}\mathit{var}(a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)}) < \infty$,

$$ \frac{m_s \mathbb{P}\mathit{var}(a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)})}{t \epsilon^2} \xrightarrow{t \to \infty} 0, $$

therefore,

$$ \sqrt{m_s t} a(J_{N_t+1}, t - S_{N_t}, V_x^{(N_t-1)}, V_x^{(N_t)}) \xrightarrow{t \to \infty} 0. $$

By Slutsky's theorem, we conclude that,

$$ \sqrt{m_s t} A_t \xrightarrow{t \to \infty} \mathcal{N}(0, \sigma_a^2). $$

To calculate $\sigma_a^2$ specifically, we have from CLT for Markov chains [52]:

$$ \begin{aligned} \sigma_a^2 &= \bar{\pi}_1 g^2 - \bar{\pi}(\bar{P}g)^2 = (\bar{\pi}_1 \bar{\pi}_2) \left( \frac{g_1^2}{g_2^2} \right) - (\bar{\pi}_1 \bar{\pi}_2) \left( \frac{(\bar{P}_{12}g_2)^2}{(\bar{P}_{21}g_1)^2} \right) \\ &= \bar{\pi}_1 g_1^2 + \bar{\pi}_2 g_2^2 - \bar{\pi}(\bar{P}_{12}g_2)^2 - \bar{\pi}_2(\bar{P}_{21}g_1)^2. \end{aligned} $$

By semi-Markov theory in (e.g., [46]), we have the following Poisson equation:

$$ g = \bar{a} + \bar{P}g, $$

with

$$ \bar{a} = a - \bar{\pi}a = a. $$

By the proof in Theorem 2, we have $\bar{\pi}a = 0$. Therefore

$$ g = a + \bar{P}g. $$

(E.1)

The solution of this equation is,

$$ g = \sum_{n \ge 0} \bar{P}^n a. $$
In matrix form

$$
\begin{align*}
g &= \begin{pmatrix} g_1 \\ g_2 \end{pmatrix} \\
&= \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} + \begin{pmatrix} 0 & \tilde{P}_{12} \\ \tilde{P}_{21} & 0 \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} + \begin{pmatrix} 0 & \tilde{P}_{12} \\ \tilde{P}_{21} & 0 \end{pmatrix}^2 \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} \begin{pmatrix} 0 & \tilde{P}_{12} \\ \tilde{P}_{21} & 0 \end{pmatrix}^3 \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} + \dots \\
&= \begin{pmatrix} a_1 + \tilde{P}_{12}a_2 + \tilde{P}_{12}\tilde{P}_{21}a_1 + \tilde{P}_{12}\tilde{P}_{21}\tilde{P}_{12}a_2 + \dots \\ a_2 + \tilde{P}_{21}a_1 + \tilde{P}_{21}\tilde{P}_{12}a_2 + \tilde{P}_{21}\tilde{P}_{12}\tilde{P}_{21}a_1 + \dots \end{pmatrix} \\
&= \begin{pmatrix} a_1 + \sum_{n>0, n \text{ odd}} \tilde{P}^n a_2 + \sum_{n>0, n \text{ even}} \tilde{P}^n a_1 \\ a_2 + \sum_{n>0, n \text{ odd}} \tilde{P}^n a_1 + \sum_{n>0, n \text{ even}} \tilde{P}^n a_2 \end{pmatrix}.
\end{align*}
$$

Therefore,

$$
\bar{P}_g = \sum_{n \ge 1} \bar{P}^n a.
$$

(a) Under high temperature

To arrive at state 1,

$$
\begin{align*}
\tilde{P}_{21} a_1 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(2, t, u, v; 1, dt', du', dv') a_1(t', u', v') \\
&= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} F_T(dt') \mathbb{1}_v(du') f_N\left(\frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} u't' = m_2 v. \\
\tilde{P}_{12} \tilde{P}_{21} a_1 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(1, t, u, v; 2, dt', du', dv') (\tilde{P}_{21} a_1)(t', u', v') \\
&= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} G(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') m_2 v' = m_2 v. \\
\tilde{P}_{21} \tilde{P}_{12} \tilde{P}_{21} a_1 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(2, t, u, v; 1, dt', du', dv') (\tilde{P}_{12} \tilde{P}_{21} a_1)(t', u', v') \\
&= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} F_T(dt) \mathbb{1}_v(du') f_N\left(\frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} m_2 v' = m_2(1 - \alpha_{\text{tan}})v. \\
\tilde{P}_{12} \tilde{P}_{21} \tilde{P}_{12} \tilde{P}_{21} a_1 &= \tilde{P}(1, t, u, v; 2, dt', du', dv') (\tilde{P}_{21} \tilde{P}_{12} \tilde{P}_{21} a_1)(t', u', v') \\
&= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} G(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') m_2 (1 - \alpha_{\text{tan}}) v' = m_2 (1 - \alpha_{\text{tan}}) v.
\end{align*}
$$

By iteration, we can conclude that,

$$
\tilde{P}^n a_1 = \begin{cases} m_2(1 - \alpha_{\text{tan}})^{\frac{n-1}{2}} v, & n \text{ odd,} \\ m_2(1 - \alpha_{\text{tan}})^{\frac{n-2}{2}} v, & n \text{ even.} \end{cases}
$$

Therefore,

$$
\sum_{n \ge 1} \tilde{P}^n a_1 = \frac{1}{\alpha_{\text{tan}}} m_2 v, \quad \text{for } n \text{ odd and even.}
$$

To arrive at state 2,

$$
\begin{align*}
\tilde{P}_{12} a_2 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(1, t, u, v; 2, dt', du', dv') a_2(t', u', v') \\
&= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} G(dt') \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') \frac{u' + v'}{2} t' = m_1 v. \\
\tilde{P}_{21} \tilde{P}_{12} a_2 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(2, t, u, v; 1, dt', du', dv') (\tilde{P}_{12} a_2)(t', u', v') \\
&= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} F_T(dt') \mathbb{1}_v(du') f_N\left(\frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} m_1 v' = m_1(1 - \alpha_{\text{tan}})v. \\
\tilde{P}_{12} \tilde{P}_{21} \tilde{P}_{12} a_2 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(1, t, u, v; 2, dt', du', dv') (\tilde{P}_{21} \tilde{P}_{12} a_2)(t', u', v') \\
&= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} G(dt') \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') m_1 v' = m_1(1 - \alpha_{\text{tan}})v. \\
\tilde{P}_{21} \tilde{P}_{12} \tilde{P}_{21} \tilde{P}_{12} a_2 &= \tilde{P}(2, t, u, v; 1, dt', du', dv') (\tilde{P}_{12} \tilde{P}_{21} \tilde{P}_{12} a_2)(t', u', v') \\
&= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} F_T(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') m_2 (1 - \alpha_{\text{tan}}) v' = m_2(1 - \alpha_{\text{tan}})^2 v.
\end{align*}
$$
By iteration, we can conclude that,

$$
\tilde{P}^n a_2 = \begin{cases} m_1(1 - \alpha_{\text{tan}})^{\frac{n-1}{2}} v, & n \text{ odd,} \\ m_1(1 - \alpha_{\text{tan}})^{\frac{n}{2}} v, & n \text{ even.} \end{cases}
$$

Therefore,

$$
\sum_{n \ge 1} \tilde{P}^n a_2 = \begin{cases} \frac{1}{\alpha_{\text{tan}}} m_1 v, & n \text{ odd,} \\ \frac{1 - \alpha_{\text{tan}}}{\alpha_{\text{tan}}} m_1 v, & n \text{ even.} \end{cases}
$$

In total,

$$
g_1 = ut + \frac{1}{\alpha_{\text{tan}}} m_1 v + \frac{1}{\alpha_{\text{tan}}} m_2 v,
$$

$$
g_2 = \frac{u+v}{2}t + \frac{1-\alpha_{\text{tan}}}{\alpha_{\text{tan}}}m_1 v + \frac{1}{\alpha_{\text{tan}}}m_2 v.
$$

So that,

$$
\tilde{P}_{21} g_1 = \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} F_T(dt') \mathbb{1}_v(du') f_N\left(\frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} \\
\times (u't' + \frac{1}{\alpha_{\text{tan}}}m_1 v' + \frac{1}{\alpha_{\text{tan}}}m_2 v')
$$

$$
= m_2 v + \frac{1 - \alpha_{\text{tan}}}{\alpha_{\text{tan}}} m_1 v + \frac{1 - \alpha_{\text{tan}}}{\alpha_{\text{tan}}} m_2 v = \frac{(1 - \alpha_{\text{tan}})m_1 + m_2}{\alpha_{\text{tan}}} v.
$$

$$
\tilde{P}_{12} g_2 = \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} G(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') \left( \frac{u' + v'}{2} t' + \frac{1 - \alpha_{\text{tan}}}{\alpha_{\text{tan}}} m_1 v' + \frac{1}{\alpha_{\text{tan}}} m_2 v' \right)
$$

$$
= m_1 v + \frac{1 - \alpha_{\text{tan}}}{\alpha_{\text{tan}}} m_1 v + \frac{1}{\alpha_{\text{tan}}} m_2 v = \frac{m_1 + m_2}{\alpha_{\text{tan}}} v.
$$

Then we are able to calculate $\sigma_a^2$,

$$
\begin{align*}
\sigma_a^2 &= \bar{\pi}_1 g_1^2 + \bar{\pi}_2 g_2^2 - \bar{\pi}(\tilde{P}_{12} g_2)^2 - \bar{\pi}_2(\tilde{P}_{21} g_1)^2 \\
&= \frac{1}{2} \mathbb{E}(X_1^2) \theta^2 + \frac{2 - \alpha_{\text{tan}}}{4} \mathbb{E}(X_2^2) \theta^2 + \frac{(m_1 + m_2)m_1}{\alpha_{\text{tan}}} \theta^2 + \frac{((1 - \alpha_{\text{tan}})m_1 + m_2)(2 - \alpha_{\text{tan}})m_2}{2\alpha_{\text{tan}}} \theta^2 \\
&= \frac{1}{2} \mathbb{E}(X_1^2) \theta^2 + \frac{2 - \alpha_{\text{tan}}}{4} \mathbb{E}(X_2^2) \theta^2 + \frac{2m_1^2 + (\alpha_{\text{tan}}^2 - 3\alpha_{\text{tan}} + 4)m_1m_2 + (2 - \alpha_{\text{tan}})m_2^2}{2\alpha_{\text{tan}}} \theta^2.
\end{align*}
$$

(b) Under low temperature

To arrive at state 1,

$$
\begin{align*}
\tilde{P}_{21} a_1 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(2, t, u, v; 1, dt', du', dv') a_1(t', u', v') \\
&= \int_{\mathbb{R} \times \mathbb{R}} \int_0^{t_0} F_T(dt') \mathbb{1}_v(du') f_N\left(\frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} u' t' \\
&\quad + \int_{\mathbb{R} \times \mathbb{R}} \int_{t_0}^{\infty} F_T(dt') \mathbb{1}_v(du') f_N(v') dv' u' t' = m_2 v.
\end{align*}
$$

$$
\tilde{P}_{12} \tilde{P}_{21} a_1 = \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(1, t, u, v; 2, dt', du', dv') (\tilde{P}_{21} a_1)(t', u', v')
$$

$$
= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} G(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') m_2 v' = m_2 v.
$$

$$
\tilde{P}_{21} \tilde{P}_{12} \tilde{P}_{21} a_1 = \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \tilde{P}(2, t, u, v; 1, dt', du', dv') (\tilde{P}_{12} \tilde{P}_{21} a_1)(t', u', v')
$$

$$
= \int_{\mathbb{R} \times \mathbb{R}} \int_0^{t_0} F_T(dt) \mathbb{1}_v(du') f_N\left(\frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} m_2 v'
$$

$$
\int_{\mathbb{R} \times \mathbb{R}} \int_{t_0}^{\infty} F_T(dt) \mathbb{1}_v(du') f_N(v') dv' m_2 v' = p(1 - \alpha_{\text{tan}}) m_2 v.
$$

$$
\tilde{P}_{12} \tilde{P}_{21} \tilde{P}_{21} a_1 = \tilde{P}(1, t, u, v; 2, dt', du', dv') (\tilde{P}_{21} \tilde{P}_{12} \tilde{P}_{21} a_1)(t', u', v')
$$

$$
= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} F(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') p(1 - \alpha_{\text{tan}}) m_2 v = p(1 - \alpha_{\text{tan}}) m_2 v.
$$
By iteration, we can conclude that,

$$ \bar{P}^n a_1 = \begin{cases} (p(1 - \alpha_{\text{tan}}))^{\frac{n-1}{2}} m_2 v, & \text{n odd,} \\ (p(1 - \alpha_{\text{tan}}))^{\frac{n-2}{2}} m_2 v, & \text{n even.} \end{cases} $$

Therefore,

$$ \sum_{n \ge 1} \bar{P}^n a_1 = \frac{1}{1 - p(1 - \alpha_{\text{tan}})} m_2 v, \quad \text{for n odd and even.} $$

To arrive at state 2,

$$ \begin{aligned} \bar{P}_{12} a_2 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \bar{P}(1, t, u, v; 2, dt', du', dv') a_2(t', u', v') \\ &= \int_{\mathbb{R} \times \mathbb{R}} \int_0^{t_0} G(dt') \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') \frac{u' + v'}{2} t' \\ &+ \int_{\mathbb{R} \times \mathbb{R} \times \mathbb{R}} \int_{t_0}^{\infty} G(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') \left( \frac{u' + v'}{2} t_0 + Y \right) f_Y(Y) dY \\ &= \left( \mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0)) t_0 \right) v. \end{aligned} $$

$$ \begin{aligned} \bar{P}_{21} \bar{P}_{12} a_2 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \bar{P}(2, t, u, v; 1, dt', du', dv') (\bar{P}_{12} a_2)(t', u', v') \\ &= \int_{\mathbb{R} \times \mathbb{R}} \int_0^{t_0} F_T(dt') \mathbb{1}_v(du') f_N\left(\frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} \\ &\quad \times \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right)v' \\ &+ \int_{\mathbb{R} \times \mathbb{R}} \int_{t_0}^{\infty} F_T(dt') \mathbb{1}_v(du') f_N(v') dv' \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right)v' \\ &= p(1 - \alpha_{\text{tan}}) \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right)v. \end{aligned} $$

$$ \begin{aligned} \bar{P}_{12} \bar{P}_{21} \bar{P}_{12} a_2 &= \int_{\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}} \bar{P}(1, t, u, v; 2, dt', du', dv') (\bar{P}_{21} \bar{P}_{12} a_2)(t', u', v') \\ &= \int_{\mathbb{R} \times \mathbb{R}} \int_0^{t_0} G(dt') \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') p(1 - \alpha_{\text{tan}}) \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right) v' \\ &+ \int_{\mathbb{R} \times \mathbb{R}} \int_{t_0}^{\infty} G(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') p(1 - \alpha_{\text{tan}}) \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right) v' \\ &= p(1 - \alpha_{\text{tan}}) \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right) v. \end{aligned} $$

$$ \begin{aligned} \bar{P}_{21} \bar{P}_{12} \bar{P}_{21} \bar{P}_{12} a_2 &= \bar{P}(2, t, u, v; 1, dt', du', dv') (\bar{P}_{21} \bar{P}_{12} a_2)(t', u', v') \\ &= \int_{\mathbb{R} \times \mathbb{R}} \int_0^{t_0} F_T(dt') \mathbb{1}_v(du') f_N\left(\frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}}\right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} \\ &\quad \times p(1 - \alpha_{\text{tan}}) \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right) v' \\ &+ \int_{\mathbb{R} \times \mathbb{R}} \int_{t_0}^{\infty} F_T(dt') \mathbb{1}_v(du') f_N(v') dv' p(1 - \alpha_{\text{tan}}) \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right) v' \\ &= (p(1 - \alpha_{\text{tan}}))^2 \left(\mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0\right) v. \end{aligned} $$

By iteration, we can conclude that

$$ \bar{P}^n a_2 = \begin{cases} (p(1 - \alpha_{\text{tan}}))^{\frac{n-1}{2}} \left( \mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0 \right) v, & \text{n odd,} \\ (p(1 - \alpha_{\text{tan}}))^{\frac{n}{2}} \left( \mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0 \right) v, & \text{n even.} \end{cases} $$

Therefore,

$$ \sum_{n \ge 1} \bar{P}^n a_2 = \begin{cases} \frac{1}{1 - p(1 - \alpha_{\text{tan}})} \left( \mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0 \right) v, & \text{n odd,} \\ \frac{p(1 - \alpha_{\text{tan}})}{1 - p(1 - \alpha_{\text{tan}})} \left( \mathbb{E}(X_1 | X_1 < t_0) + (1 - G(t_0))t_0 \right) v, & \text{n even.} \end{cases} $$
In total,

$$g_1 = ut + \frac{\mathbb{E}(X_1|X_1 < t_0) + (1 - G(t_0))t_0 + m_2}{1 - p(1 - \alpha_{\text{tan}})} v,$$

$$g_2 = \begin{cases} \frac{u+v}{2} t + \frac{p(1-\alpha_{\text{tan}})(\mathbb{E}(X_1|X_1 < t_0) + (1-G(t_0))t_0) + m_2}{1-p(1-\alpha_{\text{tan}})} v, & t < t_0 \\ \frac{u+v}{2} t_0 + Y + \frac{p(1-\alpha_{\text{tan}})(\mathbb{E}(X_1|X_1 < t_0) + (1-G(t_0))t_0) + m_2}{1-p(1-\alpha_{\text{tan}})} v, & t \ge t_0 \end{cases}$$

So that,

$$\begin{aligned} \tilde{P}_{21} g_1 &= \int_{\mathbb{R} \times \mathbb{R}} \int_0^{t_0} F_T(dt') \mathbb{1}_v(du') f_N \left( \frac{v' - (1 - \alpha_{\text{tan}})v}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} \right) \frac{dv'}{\sqrt{\alpha_{\text{tan}}(2 - \alpha_{\text{tan}})}} \\ &\quad \times (u' t' + \frac{\mathbb{E}(X_1|X_1 < t_0) + (1 - G(t_0))t_0 + m_2}{1 - p(1 - \alpha_{\text{tan}})} v') \\ &\quad + \int_{\mathbb{R} \times \mathbb{R}} \int_0^{\infty} F_T(dt') \mathbb{1}_v(du') f_N(v') dv' (u' t' + \frac{\mathbb{E}(X_1|X_1 < t_0) + (1 - G(t_0))t_0 + m_2}{1 - p(1 - \alpha_{\text{tan}})} v') \\ &= \frac{m_2 + p(1 - \alpha_{\text{tan}})(\mathbb{E}(X_1|X_1 < t_0) + (1 - G(t_0))t_0)}{1 - p(1 - \alpha_{\text{tan}})} v. \end{aligned}$$

$$\begin{aligned} \tilde{P}_{12} g_2 &= \int_{\mathbb{R} \times \mathbb{R}} \int_{t_0}^{\infty} G(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') \\ &\quad \times \left( \frac{u' + v'}{2} t_0 + Y + \frac{p(1 - \alpha_{\text{tan}})(\mathbb{E}(X_1|X_1 < t_0) + (1 - G(t_0))t_0) + m_2}{1 - p(1 - \alpha_{\text{tan}})} v' \right) f_Y(y) dy \\ &\quad + \int_{\mathbb{R} \times \mathbb{R}} \int_0^{t_0} G(dt) \mathbb{1}_v(du') \mathbb{1}_{u'}(dv') \\ &\quad \times \left( \frac{u' + v'}{2} t' + \frac{p(1 - \alpha_{\text{tan}})(\mathbb{E}(X_1|X_1 < t_0) + (1 - G(t_0))t_0) + m_2}{1 - p(1 - \alpha_{\text{tan}})} v' \right) \\ &= \frac{(\mathbb{E}(X_1|X_1 < t_0) + (1 - G(t_0))t_0) + m_2}{1 - p(1 - \alpha_{\text{tan}})} v. \end{aligned}$$

Then we are able to calculate $\sigma_a^2$, to simply the presentation, we use in the following equation,

$$m_D = \mathbb{E}(X_1|X_1 < t_0) + (1 - G(t_0))t_0.$$

Therefore,

$$\begin{aligned} \sigma_a^2 &= \tilde{\pi}_1 g_1^2 + \tilde{\pi}_2 g_2^2 - \tilde{\pi}_1 (\tilde{P}_{12} g_2)^2 - \tilde{\pi}_2 (\tilde{P}_{21} g_1)^2 \\ &= \frac{1}{2} \mathbb{E}(X_1^2) \theta^2 + \frac{m_1(m_2 + m_D)}{1 - p(1 - \alpha_{\text{tan}})} \theta^2 + \frac{1}{2} \left( \frac{m_2 + m_D}{1 - p(1 - \alpha_{\text{tan}})} \right)^2 \theta^2 + \frac{2 - \alpha_{\text{tan}}}{4} \mathbb{E}(X_2^2 | X_2 < t_0) \theta^2 \\ &\quad + \frac{1}{2} \mathbb{E}(X_2 | X_2 < t_0) \frac{p(1 - \alpha_{\text{tan}})m_D + m_2}{1 - p(1 - \alpha_{\text{tan}})} \theta^2 + \frac{1}{2} \left( \frac{p(1 - \alpha_{\text{tan}})m_D + m_2}{1 - p(1 - \alpha_{\text{tan}})} \right)^2 \theta^2 \\ &\quad + \frac{1}{4} q t_0^2 \theta^2 + \frac{q D_s}{\lambda} + \frac{1}{2} q t_0 \theta^2 - \frac{1}{2} \left( \frac{m_2 + m_D}{1 - p(1 - \alpha_{\text{tan}})} \right)^2 \theta^2 - \frac{1}{2} \left( \frac{m_2 + p(1 - \alpha_{\text{tan}})m_D}{1 - p(1 - \alpha_{\text{tan}})} \right)^2 \theta^2 \\ &= \frac{1}{2} \mathbb{E}(X_1^2) \theta^2 + \frac{m_1(m_2 + m_D)}{1 - p(1 - \alpha_{\text{tan}})} \theta^2 + \frac{2 - \alpha_{\text{tan}}}{4} \mathbb{E}(X_2^2 | X_2 < t_0) \theta^2 \\ &\quad + \frac{1}{2} \mathbb{E}(X_2 | X_2 < t_0) \frac{p(1 - \alpha_{\text{tan}})m_D + m_2}{1 - p(1 - \alpha_{\text{tan}})} \theta^2 + \frac{1}{4} q t_0^2 \theta^2 + \frac{q D_s}{\lambda} + \frac{1}{2} q t_0 \theta^2. \end{aligned}$$

In particular,

$$\mathbb{E}(X_1|X_1 < t_0) = \int_0^{t_0} t \frac{h^2}{\theta^2 t^2} e^{-\frac{h^2}{2\theta^2 t^2}} dt = \frac{h}{\theta} \sqrt{\frac{\pi}{2}} (1 - \operatorname{erf}(\frac{h}{\sqrt{2\theta t_0}})),$$

$$m_D = \frac{h}{\theta} \sqrt{\frac{\pi}{2}} (1 - \operatorname{erf}(\frac{h}{\sqrt{2\theta t_0}})) + (1 - e^{-\frac{h^2}{2\theta^2 t_0^2}}) t_0,$$

$$\mathbb{E}(X_2|X_2 < t_0) = \frac{p}{e^{-(t_0-m)/\bar{\beta}}-s} \int_0^{t_0} t dF_0(t),$$

$$\mathbb{E}(X_2^2|X_2 < t_0) = \frac{p}{e^{-(t_0-m)/\bar{\beta}}-s} \int_0^{t_0} t^2 dF_0(t).$$
## Appendix F. Proof of Theorem 4

**Proof.** If function $g: \mathbb{R} \to \mathbb{R}$ is continuous, and

$$X_n \xrightarrow[D]{t \to \infty} X,$$

then

$$g(X_n) \xrightarrow[D]{t \to \infty} g(X).$$

With CLT in Theorem 3,

$$\sqrt{m_s t A_t} = \sqrt{m_s} \frac{\Delta(X)^{(n)}}{\sqrt{t}} \xrightarrow{D} \mathcal{N}(0, \sigma_a^2).$$

We put $\Delta(X)^{(n)} = \sum_{i=1}^{n} a(J_i, X_i, V_X^{(j-1)}, V_X^{(j)})$ for the final position of particle, then we have:

$$\frac{\Delta(X)^{(n)}}{\sqrt{t}} \xrightarrow{D} \frac{\sigma_a}{\sqrt{m_s}} \mathcal{N}(0, 1).$$

If we take $g(x) = \frac{1}{2}x^2$, we have:

$$\frac{1}{2} \left( \frac{\Delta(X)^{(n)}}{\sqrt{t}} \right)^2 \xrightarrow{D} \frac{\sigma_a^2}{2m_s} \chi_1^2,$$

The diffusion coefficient of channel is then:

$$D_{\text{channel}} = \lim_{t \to \infty} \frac{(\Delta(X)^{(n)})^2}{2t} = \mathbb{E}\left(\frac{\sigma_a^2}{2m_s} \chi_1^2\right) = \frac{\sigma_a^2}{2m_s}. \quad \square$$

## Appendix G. Alternative time threshold detection and parameter estimation methods

*Generalized likelihood ratio method.* Consider that the data before the unknown time threshold follows a Fréchet distribution and after it, the distribution switches to an exponential distribution. Ignoring the distributions truncation, according to [53], we can use the generalized likelihood ratio method to detect the time threshold as follows,

$$T_{\text{detect}} = \inf \{ n \in \mathbb{N}, S_{c,n} \geq h \},$$

where we consider the comparison function:

$$S_{c,n} = \max_{1 \leq k < n} \ln \left( \frac{\sup_{\theta_1} \prod_{i=k}^{n} f_1(t_i)}{\sup_{\theta_0} \prod_{i=k}^{n} f_0(t_i)} \right).$$

Here $t_1, \dots, t_n$ are the resident times data, $\theta_j$ is the vector of all parameters under density function $f_j, j = 0, 1$. The pdf $f_0$ refers to Fréchet distribution and $f_1$ for the exponential probability density respectively. The constant $h$ is the threshold chosen, in order to satisfy $\mathbb{P}_0(T_{\text{detect}} < \infty) < \gamma$ for $\gamma \to 0$. $\mathbb{P}_0$ represents the probability distribution under the null hypothesis, which assumes that the data follows the Fréchet distribution throughout.

In this method the parameters are estimated by maximizing the likelihood. The most difficult task is to find the appropriate threshold $h$. This step is often elaborated numerically. To reduce the calculation time, as it is explained in [54], one can apply the methodology on a sliding window as follows:

$$S_{c,n} = \max_{n-m \leq k < n} \ln \left( \frac{\sup_{\theta_1} \prod_{i=k}^{n} f_1(t_i)}{\sup_{\theta_0} \prod_{i=k}^{n} f_0(t_i)} \right),$$

$$\text{such that, } \sup_{t \geq 1} \mathbb{P}_0(t \leq T_{\text{detect}} < t+m) \leq 2m/h.$$

This is a well-known and fundamental method for threshold detection. However, it remains computationally intensive when applied to large datasets.

*Q-Q plot.* Since the first part of the available data are realizations of a Fréchet distribution, a Fréchet Q-Q plot can be used to detect the time threshold. Considering the whole data or data on a sliding window $m$, to derive a Q-Q plot, parameters can be estimated via one of the following methods:

* Maximization of the log likelihood function $\sum_{i=k}^{k+m} \ln f_0(t_i)$.
<figure><img src="image_18.png" /><figcaption>Fig. H.17. Q-Q plot with different mean MC data by MLE (left) and by minimization of cdf (right).</figcaption></figure>

* Minimization of

$$e_{k,m} = \sum_{i=k}^{k+m} (\hat{F}_0(t_i) - F_0(t_i))^2, \quad 1 \le k \le N.$$

where $F_0$ is the cdf of a Fréchet distribution and $\hat{F}_0(t) = \frac{1}{N} \sum_{i=1}^{N} 1(t_i \le t)$ is the empirical cumulative distribution.

The Q-Q plot corresponding to these parameters is plotted and as soon as a deviation from the quantile plot is noticed, it is supposed that the model has switched to an exponential distribution.

**Detection by parameter estimation.** In this method, parameters of the Fréchet distribution are estimated in a sliding window of size $m$. As long as data follows a Fréchet distribution, the estimated values should not vary much with the sliding window. As soon as the estimated values switch to a different value, we can consider that the data are no longer realizations of the initial Fréchet distribution. In this framework, two estimation methods are considered.

* Maximum likelihood on a sliding window of size $m$.

* A minimization of $e_{k,m} = \sum_{i=k}^{k+m} (\hat{F}_0(t_i) - F_0(t_i))^2, \quad 1 \le k \le N$.

# Appendix H. Verification and comparison of estimation method

## H.1. Verification of estimation method

**Generalized likelihood ratio method.** For the generalized likelihood ratio method, although it is the fundamental method for transition detection, the computational time is very large. This latter makes the detection in a short delay very difficult. Moreover, due to its computational complexity, it is not the most adequate method for sensitivity analysis and asymptotic studies. That is why for the real MD data this method is not used.

**QQ-plot detection method.** The QQ-plot detection method is depicted in Figs. H.17. On the left hand side the MLE method is used to estimate parameters, and on the right hand side the cdf distance minimization is used for parameter estimation. Change point detection results in time threshold is 20 ps.

Now we explore the case that the Fréchet distribution and the exponential distribution have very close mean values. Q-Q plot is depicted in Fig. H.18. In this case the detection also gives very good results.

**Detection by parameter estimation.** Here we chose different size of windows $m \in \{1000, 3000, 5000, 10000\}$. In each window for Fréchet distribution parameter estimation, we also use two methods: window limited MLE, and minimization of window limited cdf distance. As all parameters have the same tendency in result, in this paragraph, only the behavior of $s$ estimated for different window size is depicted. In Fig. H.19 the evolution of the maximum likelihood estimate of $s$ is depicted on the left hand side. The estimates of $s$ obtained with the minimization of $e_{1,m}$ are exposed in Fig. H.19 on the right hand side. The time threshold can be evaluated to at approximately 20 ps. When the window size $m$ is too large ($m = 10000$) the results are not concluding. (See Fig. H.20.)

Afterward we explore the case that Fréchet and exponential distributions have very close mean values. The estimation $s$ by MLE is given in Fig. H.20 on the left hand side and by minimization of $e_{1,m}$ on the right hand side. It can be noticed that an abrupt large jump on parameter estimation is visible approximately at 20 ps. We can conclude that the time threshold $t_0 < 45$ ps.
<figure><img src="image_19.png" /><figcaption>Fig. H.18. Q-Q plot with similar mean MC data by MLE (left) and by minimization of cdf (right).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. H.19. Parameter *s* estimate with different mean MC data evolution by using MLE with *m* ∈ {1000, 3000, 5000, 10000} (left), Parameter *s* estimate evolution by using cdf criteria with the same window size (right).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. H.20. Parameter *s* estimate with similar mean MC data evolution by using MLE with *m* ∈ {1000, 3000, 5000, 10000} (left), Parameter *s* estimate evolution by using cdf criteria with the same window size (right).</figcaption></figure>
<table><caption>**Table H.14**
Computational time and robustness comparison of different estimation methods.</caption><thead><tr><th>Estimation data</th><th>Measurement item</th><th>Minimization of cdf</th><th>Detection by parameter estimation</th><th>QQ-plot</th></tr></thead><tbody><tr><td rowspan="2">Different mean</td><td>Mean transition time (ps)</td><td>19.9968</td><td>23.5915</td><td>19.9523</td></tr><tr><td>Transition time variance (ps)</td><td>0.0048</td><td>4.6054</td><td>0.0242</td></tr><tr><td rowspan="3">Similar mean</td><td>Mean transition time (ps)</td><td>19.9523</td><td>30.5343</td><td>19.8887</td></tr><tr><td>Transition time variance (ps)</td><td>0.0269</td><td>44.7246</td><td>0.0279</td></tr><tr><td>Computational time (s)</td><td>211.1094</td><td>2.0156</td><td>1.3594</td></tr></tbody></table>

## H.2. Comparison of estimation methods

Here, we compare the results, computational efficiency and their robustness in the follow Table H.14. The computational efficiency is measured by the time required to complete a loop of estimation. For robustness, we assess the mean and variance of the estimated transition times.

As previously mentioned, the actual transition time from MC simulation is 20 ps. For minimization of cdf, we evaluate 100 candidate values over the interval [15, 25] ps. For the parameter estimation detection method, we apply a sliding window with width $m = 1000$. Robustness is assessed by repeating each method $N = 100$ times and calculating the mean and variance of the resulting estimates.

In all cases, the transition time estimation remains consistently around 20 ps. Detection by parameter estimation gives vague graphical values rapidly, whereas QQ-plot method gives better and more robust estimation. The minimization of cdf gives a precise value, however, the interval is refined, more computational time is required. In practice, a more efficient strategy is to first identify a rough interval using QQ-plots, and then refine this interval by minimization of cdf to obtain a reliable estimation.

## Data availability

Data will be made available on request.

## References

[1] Etienne Rivard, Michel Trudeau, Karim Zaghib, Hydrogen storage for mobility: A review, Materials 12 (12) (2019) http://dx.doi.org/10.3390/ma12121973, URL https://www.mdpi.com/1996-1944/12/12/1973.

[2] Michelle Åhlén, Ocean Cheung, Chao Xu, Low-concentration CO 2 capture using metal-organic frameworks-current status and future perspectives, Dalton Trans. 52 (7) (2023) 1841–1856.

[3] Hanqian Pan, Cong Yu, Xian Suo, Lifeng Yang, Xili Cui, Huabin Xing, Emerging porous materials for carbon dioxide adsorptive capture: Progress and challenges, Mater. Chem. Front. (2023).

[4] Xiaowei Wang, Tianyi Su, Wenqing Zhang, Zhijun Zhang, Shiwei Zhang, Knudsen pumps: A review, Microsyst. Nanoeng. 6 (1) (2020) 26.

[5] Zhengyang Cai, Bilu Liu, Xiaolong Zou, Hui-Ming Cheng, Chemical vapor deposition growth and applications of two-dimensional materials and their heterostructures, Chem. Rev. 118 (13) (2018) 6091–6133.

[6] Sayan Bhowmik, Ananth Govind Rajan, Chemical vapor deposition of 2D materials: A review of modeling, simulation, and machine learning studies, Iscience 25 (3) (2022).

[7] Hongjian Wang, Meidi Wang, Xu Liang, Jinqiu Yuan, Hao Yang, Shaoyu Wang, Yanxiong Ren, Hong Wu, Fusheng Pan, Zhongyi Jiang, Organic molecular sieve membranes for chemical separations, Chem. Soc. Rev. 50 (9) (2021) 5468–5516.

[8] Rebecca L. Siegelman, Eugene J. Kim, Jeffrey R. Long, Porous materials for carbon dioxide separations, Nat. Mater. 20 (8) (2021) 1060–1072.

[9] Dong Sui, Manman Wu, Kaiyuan Shi, Changle Li, Junwei Lang, Yanliang Yang, Xiaoyan Zhang, Xingbin Yan, Yongsheng Chen, Recent progress of cathode materials for aqueous zinc-ion capacitors: Carbon-based materials and beyond, Carbon 185 (2021) 126–151.

[10] Sina Salari, Mickey Tam, Claire McCague, Juergen Stumper, Majid Bahrami, The ex-situ and in-situ gas diffusivities of polymer electrolyte membrane fuel cell catalyst layer and contribution of primary pores, secondary pores, ionomer and water to the total oxygen diffusion resistance, J. Power Sources 449 (2020) 227479.

[11] Anna Staerz, Udo Weimar, Nicolae Barsan, Current state of knowledge on the metal oxide based gas sensing mechanism, Sensors Actuators B: Chem. 358 (2022) 131531.

[12] Rouba Alrammouz, Jean Podlecki, Pascale Abboud, Brice Sorli, Roland Habchi, A review on flexible gas sensors: From materials to devices, Sensors Actuators A: Phys. 284 (2018) 209–231.

[13] Dong Feng, Zhangxin Chen, Keliu Wu, Jing Li, Xiaohu Dong, Yan Peng, Xinfeng Jia, Xiangfang Li, Dinghan Wang, A comprehensive review on the flow behaviour in shale gas reservoirs: Multi-scale, multi-phase, and multi-physics, Can. J. Chem. Eng. 100 (11) (2022) 3084–3122.

[14] Hui Wang, Li Chen, Zhiguo Qu, Ying Yin, Qinjun Kang, Bo Yu, Wen-Quan Tao, Modeling of multi-scale transport phenomena in shale gas production—a critical review, Appl. Energy 262 (2020) 114575.

[15] Cheng Liu, Zhiyong Liu, Yunsheng Zhang, A multi-scale framework for modelling effective gas diffusivity in dry cement paste: Combined effects of surface, Knudsen and molecular diffusion, Cem. Concr. Res. 131 (2020) 106035.

[16] Fumiaki Nakai, Takato Ishida, Gas diffusion in cement pastes: An analysis using a fluctuating diffusivity model, Constr. Build. Mater. 407 (2023) 133411.

[17] James Clerk Maxwell, VII. On stresses in rarified gases arising from inequalities of temperature, Philos. Trans. R. Soc. Lond. (170) (1879) 231–256.

[18] Carlo Cercignani, Maria Lampis, Kinetic models for gas-surface interactions, Transport Theory Statist. Phys. 1 (2) (1971) 101–114.
[19] R.G. Lord, Some extensions to the cercignani–lampis gas–surface scattering kernel, Phys. Fluids A: Fluid Dyn. 3 (4) (1991) 706–710.

[20] A. Filinov, M. Bonitz, D. Loffhagen, Microscopic modeling of gas-surface scattering. I. A combined molecular dynamics-rate equation approach, Plasma Sources Sci. Technol. 27 (6) (2018) 064003.

[21] A. Filinov, M. Bonitz, D. Loffhagen, Microscopic modeling of gas-surface scattering: II. Application to argon atom adsorption on a platinum (111) surface, Plasma Sources Sci. Technol. 27 (6) (2018) 064002.

[22] Johannes V. Barth, Transport of adsorbates at metal surfaces: From thermal migration to hot precursors, Surf. Sci. Rep. 40 (3–5) (2000) 75–149.

[23] T. Ala-Nissila, R. Ferrando, S.C. Ying, Collective and single particle diffusion on surfaces, Adv. Phys. 51 (3) (2002) 949–1078.

[24] Grazyna Antczak, Gert Ehrlich, Jump processes in surface diffusion, Surf. Sci. Rep. 62 (2) (2007) 39–61.

[25] Roderick Matheson Logan, R.E. Stickney, Simple classical model for the scattering of gas atoms from a solid surface, J. Chem. Phys. 44 (1) (1966) 195–201.

[26] Roderick Matheson Logan, James C. Keck, Classical theory for the interaction of gas atoms with solid surfaces, J. Chem. Phys. 49 (2) (1968) 860–876.

[27] Natasa Mateljevic, Jay Kerwin, Sharani Roy, J.R. Schmidt, John C. Tully, Accommodation of gases at rough surfaces, J. Phys. Chem. C 113 (6) (2009) 2360–2367.

[28] Kyoji Yamamoto, Hideki Takeuchi, Toru Hyakutake, Scattering properties and scattering kernel based on the molecular dynamics analysis of gas-wall interaction, Phys. Fluids 19 (8) (2007).

[29] Zijing Wang, Chengqian Song, Fenghua Qin, Xisheng Luo, Establishing a data-based scattering kernel model for gas-solid interaction by molecular dynamics simulation, J. Fluid Mech. 928 (2021) A34.

[30] Shahin Mohammad Nejad, Eldhose Iype, Silvia Nedea, Arjan Frijns, David Smeulders, Modeling rarefied gas-solid surface interactions for couette flow with different wall temperatures using an unsupervised machine learning technique, Phys. Rev. E 104 (1) (2021) 015309.

[31] Tengfei Liang, Qi Li, Wenjing Ye, A physical-based gas-surface interaction model for rarefied gas flow simulation, J. Comput. Phys. 352 (2018) 105–122.

[32] Tengfei Liang, Jun Zhang, Qi Li, A parameter-free physical model for gas-surface interaction, Phys. Fluids 33 (8) (2021).

[33] Pierre Magnico, Quy-Dong To, Collisions, adsorption and self diffusion of gas in nanometric channels by molecular dynamics and stochastic simulation and the case of helium gas in graphitic slit pore, Int. J. Heat Mass Transfer 214 (2023) 124371.

[34] M. v. Smoluchowski, Zur kinetischen theorie der transpiration und diffusion verdünnter gase, Ann. Phys., Lpz. 338 (16) (1910) 1559–1570.

[35] Gaurav Arya, Hsueh-Chia Chang, Edward J. Maginn, Knudsen diffusivity of a hard sphere in a rough slit pore, Phys. Rev. Lett. 91 (2) (2003) 026102.

[36] Nikolaos Limnios, Gheorghe Oprisan, Semi-Markov Processes and Reliability, Springer Science & Business Media, 2012.

[37] Andrzej Świderski, Anna Borucka, Małgorzata Grzelak, Leszek Gil, Evaluation of machinery readiness using semi-Markov processes, Appl. Sci. 10 (4) (2020) 1541.

[38] Jacques Janssen, Semi-Markov models in economics and insurance, in: F. de Vylder, M. Goovaerts, J. Haezendonck (Eds.), Premium Calculation in Insurance, Springer Netherlands, Dordrecht, 1984, pp. 241–261, http://dx.doi.org/10.1007/978-94-009-6354-2_14.

[39] Nikolaos Limnios, Anatoliy Swishchuk, Discrete-time semi-Markov random evolutions and their applications, Adv. in Appl. Probab. 45 (1) (2013) 214–240.

[40] Irene Votsi, Nikolaos Limnios, George Tsaklidis, Eleftheria Papadimitriou, Estimation of the expected number of earthquake occurrences based on semi-Markov models, Methodol. Comput. Appl. Probab. 14 (2012) 685–703.

[41] J.H.A. de Smit, The single server semi-markov queue, Stochastic Process. Appl. 22 (1) (1986) 37–50, http://dx.doi.org/10.1016/0304-4149(86)90112-2.

[42] Vlad Stefan Barbu, Nikolaos Limnios, Semi-Markov Chains and Hidden Semi-Markov Models Toward Applications: Their Use in Reliability and DNA Analysis, volume 191, Springer Science & Business Media, 2009.

[43] Bei Wu, Lirong Cui, Chen Fang, Reliability analysis of semi-Markov systems with restriction on transition times, Reliab. Eng. Syst. Saf. 190 (2019) 106516, http://dx.doi.org/10.1016/j.ress.2019.106516.

[44] Fei Liu, Shihao Xia, Shanhe Su, Semi-Markov processes in open quantum systems. III. Large deviations of first passage time statistics, 2024, arXiv preprint arXiv:2407.01940.

[45] Patrick Billingsley, Convergence of Probability Measures, John Wiley & Sons, 2013.

[46] Søren Asmussen, Applied Probability and Queues, volume 2, Springer, 2003.

[47] Haruhisa Nishino, Kazuhiko Kakamu, Grouped data estimation and testing of Gini coefficients using lognormal distributions, Sankhya B 73 (2011) 193–210.

[48] Chris Bambey Guure, Noor Akma Ibrahim, Mohd Bakri Adam, et al., Bayesian inference of the Weibull model based on interval-censored survival data, Comput. Math. Methods Med. 2013 (2013).

[49] Hatim Solayman Migdadi, Mohammad Subhi Al-Batah, Bayesian inference based on the interval grouped data from the weibull model with application, Br. J. Math. Comput. Sci. 4 (9) (2014) 1170–1183.

[50] Neil A. Weiss, Paul T. Holmes, Michael Hardy, A course in probability, Pearson Addison Wesley Boston, MA, USA, 2006.

[51] Iosif Ilich Gikhman, Anatolii Vladimirovich Skorokhod, The Theory of Stochastic Processes I, Springer-Verlag, 1974, URL https://nla.gov.au/nla.cat-vn291069.

[52] A. Gut, Stopped random walks: Limit theorems and applications, in: Springer Series in Operations Research and Financial Engineering, Springer New York, 2009, URL https://books.google.fr/books?id=twFBs5G7jHIC.

[53] Michele Basseville, Igor V. Nikiforov, Detection of Abrupt Changes: Theory and Application, volume 104, Prentice hall Englewood Cliffs, 1993.

[54] Tze Leung Lai, Information bounds and quick detection of parameter changes in stochastic systems, IEEE Trans. Inform. Theory 44 (7) (1998) 2917–2929.