# Weak dissipation drives and enhances Wada basins in three-dimensional
chaotic scattering 111This work has been supported by the Spanish State
Research Agency (AEI) and the European Regional Development Fund (ERDF, EU)
under Project No. PID2019-105554GB-I00.

Diego S. Fernández  Jesús M. Seoane<EMAIL_ADDRESS>Miguel A.F. Sanjuán
Nonlinear Dynamics, Chaos and Complex Systems Group, Departamento de Física,
Universidad Rey Juan Carlos
Tulipán s/n, 28933 Móstoles, Madrid, Spain  Department of Applied Informatics,
Kaunas University of Technology
Studentu 50-415, Kaunas LT-51368, Lithuania

###### Abstract

Chaotic scattering in three dimensions has not received as much attention as
in two dimensions so far. In this paper, we deal with a three-dimensional open
Hamiltonian system whose Wada basin boundaries become non Wada when the
critical energy value is surpassed in the absence of dissipation. In
particular, we study here the dissipation effects on this topological change,
which has no analogy in two dimensions. Hence, we find that non-Wada basins,
expected in the absence of dissipation, transform themselves into partially
Wada basins when a weak dissipation reduces the system energy below the
critical energy. We provide numerical evidence of the emergence of the Wada
points on the basin boundaries under weak dissipation. According to the paper
findings, Wada basins are typically driven, enhanced and, consequently,
structurally stable under weak dissipation in three-dimensional open
Hamiltonian systems.

###### keywords:

Chaotic scattering , Open Hamiltonian system , Three-dimensional system ,
Dissipative system , Wada property , Merging method

##  1 Introduction

Chaotic scattering [1] is an important field in nonlinear dynamics with a wide
range of different applications in physics, such as fluid dynamics [2] or
material transitions [3], but also in biology [4] or medicine [5], to name a
few. A chaotic scattering process is a direct manifestation of transient chaos
[6] in open dynamical systems. In these processes, which can be modeled by
means of open Hamiltonian systems, typically a point particle interacts with a
nonlinear potential function for a finite time, or transient, until it escapes
through one of the possible exits to infinity. The exit basins, sets of
initial conditions that lead to a certain exit, are a common tool in this kind
of systems to ascertain which initial condition regions are unpredictable,
since a minimum uncertainty in the initial conditions can hinder the exit
prediction. The vast majority of works on open Hamiltonian systems have been
devoted to conservative two-dimensional problems [7, 8]. Nonetheless, these
systems have also been studied in the presence of external perturbations, such
as periodic driving [9], additive noise [10] or weak dissipation [11, 12, 13].

Specifically, in the presence of dissipation, the unpredictability in
determining the system final state has been studied, both in discrete maps
[11] and two-dimensional open Hamiltonian systems [12, 13]. These works
conclude that the nonhyperbolic dynamics, dominated by Kolmogorov-Arnold-Moser
(KAM) tori, is suddenly destroyed by introducing an even extremely small
amount of dissipation. As a consequence of this, the more intense the
dissipation rate, the less the system unpredictability. On the other hand, the
Wada property has been also investigated in two-dimensional open Hamiltonian
systems, whose exit basins typically exhibit this property [14]. Furthermore,
Wada basins have been demonstrated to be structurally stable in the presence
of weak dissipation in these systems [12]. This is an intriguing topological
property of some sets, according to which the boundary separating three or
more sets is common to all of them [15]. In this manner, on a Wada boundary is
always possible to find all other basins on the boundary between two of them,
which represents a clear source of system unpredictability. For further
details, see Refs. [16, 17], where the state of the art of the Wada property
and the current methods to detect it are reviewed.

However, three-dimensional open Hamiltonian systems under dissipation still
require to be investigated, since the conservative scenario is the one that
has been the focus of attention so far [18, 19, 20, 21]. In this paper, we
deal with a three-dimensional open Hamiltonian system with four possible exits
[22], which has a two-dimensional version [23]. In particular, in the absence
of dissipation, a striking basin topological metamorphosis takes place with no
analogy in two dimensions: Wada basins become non-Wada basins when the system
energy exceeds a certain critical value $E_{c}$. In this context, we wonder
about the topological consequences of introducing dissipation, since the
energy loss can determine whether the basins are Wada or not.

The organization of this paper is as follows. In Sec. 2, we describe the
Hamiltonian dynamics of the model in the conservative and the dissipative
cases. In Sec. 3, we compute the system unpredictability versus the
dissipation parameter by means of the fractal dimension of the chaotic saddle.
This fractal dimension is key to understand the Wada basins in the dissipative
system. In Sec. 4, we propose the dissipative mechanism by which non-Wada
basins in the conservative case become Wada. Furthermore, we compute the
fraction of Wada points on the basin boundaries as the dissipation parameter
is varied. Finally, in the conclusions, we summarize the main findings and
relate them to possible applications where they may be of interest.

##  2 Model description

In order to study the dissipation effects on a three-dimensional chaotic
scattering problem, we use a system that models the interaction of a point
particle with non-rotating diatomic molecules [22]. The system dynamics is
modeled by the dimensionless and time-independent Hamiltonian

$H(p_{x},p_{y},p_{z},x,y,z)=\frac{p_{x}^{2}+p_{y}^{2}+p_{z}^{2}}{2}+V(x,y,z),$
(1)

where $p_{x}$, $p_{y}$ and $p_{z}$ are the momenta and $V(x,y,z)$ is a
nonlinear potential function which depends on the positions $x$, $y$ and $z$.
The potential consists of four Morse hills [24] whose centers are placed at
the vertices of a regular tetrahedron of unit side lengths. Specifically,
located at $(x_{1},y_{1},z_{1})=(0,0,\sqrt{2/3})$,
$(x_{2},y_{2},z_{2})=(1/2,-1/(2\sqrt{3}),0)$,
$(x_{3},y_{3},z_{3})=(-1/2,-1/(2\sqrt{3}),0)$ and
$(x_{4},y_{4},z_{4})=(0,1/\sqrt{3},0)$, depicted as red dots in Fig. 1. For
clarity, the tetrahedron centroid, $(x_{c},y_{c},z_{c})=(0,0,1/\sqrt{24})$, is
also depicted as a black dot in Fig. 1. Here the three-dimensional Morse hills
have spherical symmetry but slightly distorted by the presence of the others.
Each one of them consists of a repulsive core (any tetrahedron vertex)
surrounded by an attractive region, which is less effective, and ultimately
negligible, if the particle moves away from the hill center. The Morse
potential represents short-range interactions, implying that the system
chaotic behavior is concentrated near the tetrahedron formed by the four
hills’ centers. In this manner, we define the scattering region as the space
bounded by the four tetrahedron planes. That is, we consider that a particle
escapes when it crosses any of these planes. The complete potential function
is written as

$V(x,y,z)=\sum_{i=1}^{4}V_{i}(x,y,z),$ (2)

where

$V_{i}(x,y,z)=\frac{V_{0}}{2}\left[1-e^{-\alpha(r_{i}-r_{e})}\right]^{2}-\frac{V_{0}}{2}$
(3)

represents the potential function associated with each one of the four hills,
and

$r_{i}(x,y,z)=\sqrt{(x-x_{i})^{2}+(y-y_{i})^{2}+(z-z_{i})^{2}}$ (4)

represents the distances between the particle position and each one of the
four hills’ centers. In addition, the quantities $V_{0}$, $\alpha$ and $r_{e}$
are the Morse potential parameters, which are related to the well depth, the
well steepness and the effective range of each hill, respectively. Without
loss of generality, we set along this paper the parameter values as $V_{0}=1$,
$\alpha=6$ and $r_{e}=0.68$.

Given a system energy value $E$, there exists a three-dimensional region of
the physical space inaccessible to the particle, where $E<V$ is satisfied.
That is, the inaccessible region limits are simply the potential barriers
where the particle usually bounces off before escaping. In this specific
system, when the energy is greater than a certain critical energy $E_{c}$, the
inaccessible region is formed by four disconnected spherical-like regions
surrounding the Morse hills’ centers, as shown in Fig. 1(a). As the system
energy is decreased, the four disconnected regions grow in radius, touching
each other when the critical energy is reached. Finally, when the energy is
smaller than the critical energy, the inaccessible region is a single fully
connected region, as observed in Fig. 1(b). In addition, we display the
potential function in the plane $z=0$ in Figs. 1(c) and 1(d) to better
visualize how inaccessible regions remain connected or disconnected depending
on the energy. The existence of this critical energy is key to the phase space
topology, even in presence of dissipation, as we explain in detail below. We
compute $E_{c}\approx 2.250414$ by using the critical trajectory concept,
which has been developed in two-dimensional open Hamiltonian systems [25].
Hence, we work along the paper on a range of initial energy values
$E_{0}\in[2,3]$, which includes the critical energy.

Furthermore, here we focus on the dissipative case, in which the system
initial energy $E_{0}$ decreases as the time passes by. As is common in this
kind of problems, we introduce the dissipation parameter $\mu$ in Hamilton’s
equations of motion by adding a term proportional to the particle velocity in
the three degrees of freedom [26], yielding

$\displaystyle\dot{x}=p_{x},\quad\dot{p_{x}}=-\frac{\partial V}{\partial
x}-\mu p_{x},$ (5)
$\displaystyle\dot{y}=p_{y},\quad\dot{p_{y}}=-\frac{\partial V}{\partial
y}-\mu p_{y},$ $\displaystyle\dot{z}=p_{z},\quad\dot{p_{z}}=-\frac{\partial
V}{\partial z}-\mu p_{z}.$

Then, five attractors appear when the dissipation is introduced to which the
particle may converge and be trapped forever. Specifically, one of them
appears at the tetrahedron centroid, whereas the remaining four attractors are
associated with potential relative minima, depicted as blue dots in Fig. 1. As
four attractors lie outside the defined scattering region, we take into
account in the dissipative case that a particle escapes when it moves an
arbitrary distance $l$ away from the tetrahedron centroid, e.g., $l>3$, such
that it is highly unlikely that the particle can be trapped in any attractor.
On the other hand, Fig. 2(a) shows how the escape time average is affected by
the dissipation. We notice that, from $\mu>10^{-2}$, the escape time averages
in the dissipative and conservative cases begin to differentiate and from such
a value the introduced dissipation becomes relevant. So, we use the value
range $\mu\in[0.01,0.5]$ along the present work.

Three examples of trajectories are plotted in Fig. 2(b). In particular, they
evolve for the same dissipation rate but depart from different values of
$E_{0}$. Interestingly, the trajectory starting from the lowest initial energy
manages to escape from the scattering region in the first term, but is
eventually trapped in one of the attractors. Conversely, the other two
trajectories escape quickly from the tetrahedron barely interacting with the
potential. This last kind of behavior is the one we observe the most in the
computations under weak dissipation.

Finally, the conservative (i.e., $\mu=0$) and dissipative dynamics is
integrated by means of an adaptive Runge-Kutta-Fehlberg method [27] with a
relative tolerance of $10^{-12}$ for each equation system solution. The
maximum value of the integration step is $5\times 10^{-3}$, small enough to
conserve the mechanical energy along the particle trajectory in the
conservative case.

Figure 1: The physical space for (a) $E=3>E_{c}$ and (b) $E=2<E_{c}$, where
$E_{c}\approx 2.25$. The interior of the gray spherical-like regions
represents the physical space part inaccessible to the particle, which remain
(a) disconnected for $E>E_{c}$ and (b) connected for $E<E_{c}$. We define the
scattering region as the regular tetrahedron (dashed red lines) formed by the
four Morse hills’ centers (red dots), since the chaotic behavior is only
relevant near the tetrahedron centroid (black dot). Five attractors (black and
blue dots) appear when the dissipation is introduced, and the particle may
trap in them forever. Finally, we show a potential function section with the
plane $z=0$ for (c) $E=3$ and (d) $E=2$, where the minimum value is the same
in both cases.

Figure 2: (a) Escape time average as the dissipation parameter $\mu$ is varied
for fifty different values belonging to the interval $\mu\in[10^{-4},0.5]$.
The three curves correspond to three initial energy cases $E_{0}=2<E_{c}$
(blue), $E_{0}=2.5>E_{c}$ (red) and $E_{0}=3\gg E_{c}$ (green), where
$E_{c}\approx 2.25$. The horizontal dashed lines represents the escape time
average in the conservative case. The vertical dashed line is placed at
$\mu=10^{-2}$ and indicates the beginning of the regime where the dissipation
significantly affects the system dynamics. The escape times have been computed
by launching 2500 particles from the plane $z=0$ toward the tetrahedron
centroid and taking into account if $l>3$. Finally, the three black circles
correspond to the parameter values of the trajectories shown in (b). These
trajectories are launched from the same initial condition (square) and evolve
for $\mu=0.3$ and $E_{0}=2$ (blue), $E_{0}=2.5$ (red) and $E_{0}=3$ (green),
respectively. Finally, we depict the five attractors and the scattering region
(black).

##  3 Non-monotonic behavior of the chaotic saddle dimension under dissipation

Along this section, we use the exit basins to analyze the system
unpredictability. Specifically, we choose the plane $z=0$ to build the basins.
In addition, we launch every particle toward the tetrahedron centroid to
simulate trajectories that capture all the three-dimensional scattering region
phenomenology. Four basins (or colors) appear, which correspond to the four
tetrahedron planes through which the particle can escape. The yellow basin is
associated with the bottom exit, i.e., the plane $z=0$, whereas the other
basins (red, green and blue) are associated with the three remaining
tetrahedron planes. Figure 3 shows the exit basins in the conservative case
for several energy values. As expected, the higher the energy, the less
fractal the basin boundaries. The potential function exits become wider as the
energy is increased, implying the escape to happen without spending long
transients. As a matter of fact, trajectories can amplify extremely small
differences in their initial conditions during long transients due to the
nonlinear dynamics, and eventually escape through different exits. For this
reason, fractal basin boundaries are a well-known source of the
unpredictability of the system final state.

Figure 3: Exit basins in the conservative case for (a) $E=2<E_{c}$, (b)
$E=2.26\approx E_{c}$, (c) $E=2.5>E_{c}$, and (d) $E=3\gg E_{c}$, where
$E_{c}\approx 2.25$. The yellow color corresponds to particles that escape
through the bottom plane, whereas the other colors (red, blue and green) are
related to the remaining planes. Finally, the white color represents the
inaccessible initial conditions, which satisfy $E<V(x,y,0)$.

The nonattracting invariant set called chaotic saddle [28] is responsible for
the appearance of the fractal basin boundaries. The numerical computation of
the fractal dimension $D_{c}$ of the chaotic saddle is a common method to
quantify this unpredictability in chaotic scattering problems [29]. For
instance, the bigger the computed $D_{c}$ value, the more unpredictable the
system, since a bigger region the phase space (the exit basins) is occupied by
the chaotic saddle (the fractal basin boundaries). In the A at the end of the
manuscript we derive the expression

$D_{c}=2d+1,$ (6)

which relates $D_{c}$ to the fractal dimension $d$ of the exit basin
boundaries. Thus, we can determine $D_{c}$ after calculating $d$ with the help
of the uncertainty exponent algorithm [30] (see the A for further computation
details).

Therefore, we display the fractal dimension of the chaotic saddle versus the
initial energy and the dissipation in Figs. 4(a) and 4(b), respectively. On
the one hand, the greater the initial energy, the smaller the value of
$D_{c}$. As already mentioned, increasing the energy softens the basin
boundaries, making the system to be more predictable. On the other hand, the
dissipative case is more complicate. When the dissipation is weak, the fractal
dimension $D_{c}$ increases until a maximum value of $D_{c}$ is obtained, as
observed for three different initial energies in Fig. 4(b). Immediately after
reaching this maximum value, the system unpredictability exhibits an abrupt
decreasing to presumably $D_{c}=3$, as theoretically reasoned in the A. This
tendency of $D_{c}$ can be interpreted by visualizing the dissipation
parameter impact on the basin boundaries.

Figure 4: (a) Fractal dimension of the chaotic saddle as $E_{0}$ is varied for
three different dissipation values. The vertical dashed line is placed at
$E_{c}\approx 2.25$, and four circles are drawn at the energy values of the
exit basins already shown in Fig. 3. (b) Similarly, the fractal dimension of
the chaotic saddle as $\mu$ is varied for three different initial energies.
The vertical dashed lines are placed at $\mu_{m}(E_{0})$, i.e., the
dissipation parameter values where the fractal dimension $D_{c}$ reaches its
maximum for a given initial energy. Specifically, $\mu_{m}(2)=0.33$,
$\mu_{m}(2.5)=0.39$ and $\mu_{m}(3)=0.45$. Again, we draw four circles at the
$\mu$ values of the exit basins shown in Fig. 5. See the A for more
computation details.

Hence, Fig. 5 shows four dissipative cases for $E_{0}=2.5$. Firstly, the
basins are barely affected by the dissipation for $\mu=0.01$. However, the
higher the dissipation, the longer the transients within the scattering
region. Then, as observed for $\mu=0.2$, the basin boundaries widen and the
value of $D_{c}$ increases. A weak dissipation makes the system to be more
unpredictable. In addition, as the dissipation is increased, some trajectories
starting from the basin boundaries can be trapped. Then, another basin (or
color) appears and the basins become more complex and unpredictable. In this
situation, the fractal dimension $D_{c}$ grows until a maximum value is
reached when all trajectories starting from the boundaries are trapped. The
latter is satisfied for $\mu=0.39$, where the fractal basin boundaries are
completely erased by a new white and smooth basin. Here we can define
$\mu_{m}$ as the dissipation parameter value that causes the fractal dimension
$D_{c}$ to be maximum. In particular, the value of $\mu_{m}$ depends on the
initial energy, since the larger $E_{0}$ the more dissipation must be
introduced to trap particles. Finally, the new white basin grows when the
dissipation is stronger than $\mu_{m}$, as seen for $\mu=0.45$. Therefore, the
fractal dimension $D_{c}$ sharply declines in this scenario of strong
dissipation.

Figure 5: Exit basins in the dissipative case for $E_{0}=2.5$ and (a)
$\mu=0.01\ll\mu_{m}$, (b) $\mu=0.2<\mu_{m}$, (c) $\mu_{m}=0.39$, and (d)
$\mu=0.45>\mu_{m}$. Specifically, $\mu_{m}=0.39$ represents the dissipation
value for which the fractal dimension is maximum for $E_{0}=2.5$. (Inset plot)
We magnify a basin region to illustrate how the fractal basin boundary at
$\mu_{m}=0.39$ is replaced by a new white basin, which represents the initial
conditions that lead trajectories to be trapped in the attractors.

##  4 Driving and enhancement of Wada basins under dissipation

Now we investigate how dissipation affects the topology of the basin
boundaries, which may exhibit the Wada property. Briefly, Wada basins share a
common boundary implying that between two basins it is always possible to find
all other basins. The exit basin boundary is non Wada for $E>E_{c}$ in the
conservative regime. That is, only some isolated points on the boundaries are
Wada points, which are defined as points separating the four basins (or
colors). For $E>E_{c}$, the inaccessible regions remain disconnected, causing
a part of the tetrahedron edges to be accessible. The exits, i.e., the
tetrahedron planes, are connected by the part of the edges not swallowed by
the inaccessible regions. Hence, it is possible to propose a continuous and
smooth path of initial conditions whose corresponding trajectories escape
through an exit at the path beginning and through another different exit at
the path ending. Then, necessarily, a trajectory must escape at some path
point through the tetrahedron edge separating the two exits. This hypothetical
path of initial conditions consists of two one-dimensional basins separated by
a boundary where only two colors meet, i.e., a non-Wada point.

On the other hand, the exit basin boundary exhibits the Wada property for
$E<E_{c}$, i.e., all boundary points are Wada points. The inaccessible region
is connected and completely covers the tetrahedron edges. Hence, the exits
remain disconnected. In this case, there is no hypothetical path of initial
conditions whose corresponding trajectories escape through only two exits.
These trajectories, as they approach the now inaccessible edge of the
tetrahedron, deviate until escaping inevitably through a third exit different
from the previous ones. Note that this last argument is valid for any path
separating two basins, no matter its length. Therefore, a different basin
always appears between any two basins, and so on. For $E<E_{c}$, the path is
formed by two one-dimensional basins separated by a Wada boundary, consisting
of an infinite number of Wada points.

Figure 6: Energy values of fifty particles over time for three dissipative
cases, namely, $E_{0}=2.5$ and $\mu=0.15$ (blue), $E_{0}=2.5$ and $\mu=0.35$
(yellow), and $E_{0}=2.75$ and $\mu=0.15$ (green). We draw a red (black) dot
if the particle final energy before leaving the tetrahedron is smaller
(greater) than the critical energy $E_{c}\approx 2.25$, which is indicated by
the horizontal red dashed line.

In this context, where the system energy determines whether the basins are
Wada or not, introducing energy dissipation may affect their topology. For
instance, as illustrated in Fig. 6, a particle starting from $E_{0}>E_{c}$ may
leave the tetrahedron with $E<E_{c}$ due to dissipation after spending a
finite time within the scattering region. Hence, this particle escapes from
the tetrahedron with the inaccessible regions fully connected. Regarding the
latter, it is worth recalling that a dissipative system can be understood as a
sequence of momentarily conservative systems of decreasing energy over time.
So, the momentarily inaccessible regions of the potential function grow in
radius as the energy is dissipated. Then, when all particles starting from the
basin boundaries escape with energy values lower than the critical energy,
only Wada paths of initial conditions are possible. In this manner, from a
theoretical point of view, non-Wada basins expected in the conservative case
may become Wada basins under dissipation. Nonetheless, when the dissipation
allows only a few particles to escape with energies below the critical energy,
the basin boundaries are formed by Wada and non-Wada points, which is known as
the partially Wada property [31]. In this situation, only trajectories with
long transients contribute to the formation of Wada points. Furthermore, as
detailed in the previous section, a strong dissipation is expected to reduce
both the fractal boundary and the number of Wada points.

Now we test the emergence of the Wada points in the dissipative case by
applying a recently developed numerical method called merging method [32]. The
method is based on a non-obvious aspect of Wada basins: if the boundary of a
basin is Wada, then this boundary remains invariant under the merging of the
other basins. Thus, the method consists of verifying this latter statement by
merging basins. On the other hand, an advantage of this method is that it only
requires a basin image of finite resolution as input. This is useful when
there is no more knowledge about the system and its dynamics than a basin
image. However, the method has also limits, since it only verifies whether a
boundary point is Wada or non-Wada up to a given resolution. In addition, the
Wada point detection depends strongly on the fattening parameter $r$. The
method only ascertains what the fraction of Wada points is given a value of
$r$. In practice, the computation has to be repeated increasing
$r\in\mathbf{N}$ to test how the fraction of Wada points evolves as a function
of $r$. Finally, the computation ends at the value of $r$ for which the Wada
point fraction as a function of $r$ has converged. The steps and many other
aspects of this method are extensively detailed in Ref. [32].

Figure 7: (a) Exit basins as computed for $E_{0}=2.5$ and $\mu=0.01$ with a
1000$\times$1000 resolution. We focus on the magnified region
$x_{0}\in[-0.029,0)$ and $y_{0}\in[-0.249,-0.19]$, because non-Wada points can
be detected when considering these parameter values and finite resolution.
Then, here we implement the merging method to compute the Wada point fraction
as the dissipation parameter is varied [32]. (b) The four merged basins are
represented in gray and are obtained by merging the remaining basins
(including the white basin) except the colored one, namely, yellow, blue, red
and green, respectively.

We provide an example of merged basins for $E_{0}=2.5$ in Fig. 7. Here we only
focus on a convenient basin region where non-Wada points can be detected when
considering a finite resolution for $E_{0}=2.5$. As a matter of fact, the vast
majority of boundary points seem to be Wada at lower resolutions despite
$E>E_{c}$ (e.g., see again Fig. 3(c)). We have verified by successive
magnifications of the basin boundary that, at higher resolution, the boundary
points only separate two basins for $E>E_{c}$ in the conservative case, except
for some isolated Wada points where the four colors meet. This kind of fractal
structure is known as a slim fractal [33], which is simpler at smaller
resolutions.

Figure 8: (a) The fraction of Wada points $f_{W}$ on the basin boundaries for
several values of the fattening parameter and (b) the fraction of trapped
particles $f_{t}$ both versus one hundred values of the dissipation parameter
belonging to the interval $\mu\in[0.01,0.45]$ and $E_{0}=2.5$. We draw a red
dashed line at the maximum value of $f_{W}$, which is reached when dissipation
begins to be strong enough to trap particles, as indicated in (b) as well. On
the other hand, $\mu_{m}=0.39$ (black dashed line) represents strong
dissipation for $E_{0}=2.5$ and, regarding the basin boundary topology, it
indicates the separation between the partially Wada [31] and non-Wada regimes.
Furthermore, we draw four circles at the values of $\mu$ of the basins shown
in Fig. 9.

The fractions of Wada points $f_{W}$ and trapped particles $f_{t}$ are
displayed versus the dissipation parameter in Figs. 8(a) and 8(b),
respectively. In light of the results, the method converges for $r=6$.
Firstly, note that the value of $f_{W}$ grows slightly when the dissipation is
weak enough that the trapped particle fraction is negligible, e.g., roughly,
$f_{t}<10^{-3}$. Here the proposed mechanism drives the formation of new
detectable Wada points and $f_{W}$ reaches a maximum value, as indicated by a
red dashed line in both figures. However, as the dissipation is increased and
strong enough to trap a substantial fraction of particles, the number of Wada
points decreases drastically. In particular, as already explained above, a
white basin appears and the probability of detecting Wada points decreases,
since a colored boundary point is less likely to be surrounded by the other
three remaining colors. Finally, the fraction of Wada points tends to vanish,
$f_{W}\to 0$, when the dissipation exceeds $\mu_{m}=0.39$, which is considered
as strong dissipation for $E_{0}=2.5$. We indicate the beginning of the strong
dissipative regime by a black dashed line in both figures. This result is in
perfect agreement with the findings obtained in Sec. 3, where we prove for
$\mu>\mu_{m}$ that all fractal boundaries are removed and replaced by a new
smooth basin representing trapped particles. Hence, for $\mu>\mu_{m}$ all
detected boundary points are non-Wada, so that $f_{W}=0$. The whole process of
basin transformation described here is illustrated in Fig. 9, where the
appearance of Wada points under weak dissipation and their drastic destruction
in the strong dissipation regime can be appreciated.

Figure 9: Exit basins (gray scale) as computed for $E_{0}=2.5$, and (a)
$\mu=0.09$, (b) $\mu=0.24556$, (c) $\mu=0.36111$ and (d) $\mu=0.41$. We plot
in red the Wada points detected by the merging method when the fattening
parameter value is $r=6$. Finally, the white basin represents particles
trapped in the attractors.

##  5 Conclusions

Along this paper, we have studied the dissipation effects on a three-
dimensional open Hamiltonian system. Specifically, in the conservative regime,
the basin boundaries undergo a topological metamorphosis at the critical
energy $E_{c}\approx 2.25$. Then, the Wada basins become non-Wada basins when
this critical energy is surpassed. In this context, we argue that non-Wada
basins expected in the conservative case (for $E>E_{c}$) can transform
themselves into partially Wada basins in some situations where dissipation is
able to reduce the system energy below $E_{c}$. As a matter of fact, we
confirm the latter and detect how Wada points emerge on the basin boundaries
under weak dissipation, whereas strong dissipation prevents the formation of
Wada points. This phenomenology has no counterpart in dissipative two-
dimensional open Hamiltonian systems, because there exists no critical
energies related to a basin topological metamorphosis. As a result of this
research, we claim that the Wada property is driven, enhanced and,
consequently, structurally stable under weak dissipation in three-dimensional
open Hamiltonian systems.

The analysis carried out here may be useful in other three-dimensional models
in which the aim is to increase their predictability. Since a control
parameter might exceed a certain threshold value, an unpredictability source
such as a Wada boundary can be eliminated. However, we have proved that this
control parameter can also be modified by the energy loss over time. In this
situation, the appearance of Wada basins would again hinder the
predictability. For instance, an application where the friction has crucial
consequences might be the three-dimensional chaotic open flows of active
particles. They can model from the dynamics of pollutants in the atmosphere
(chemical activity) to the plankton population dynamics in an ocean
(biological activity) [34]. Moreover, there exists a rising interest in soft
Lorentz gases because they are able to simulate electronic transport in
graphene-like structures [35]. Diffusion phenomenon in such systems can be
modeled by means of three-dimensional open Hamiltonian systems arranged in
complex lattices, so a weak dissipation can hinder their predictability.
Finally, it is worth mentioning that dissipation is not the only factor that
can drive and enhance Wada structures. In fact, some systems have been proven
to exhibit the Wada property when time delays are involved [36]. Thus,
dissipation is one more factor to be taken into account when predicting the
behavior of dynamical systems, specially, in three dimensions.

##  Appendix A Derivation and computation of the fractal dimension $D_{c}$ of
the chaotic saddle.

Since the chaotic saddle is the intersection of the stable and unstable
manifolds [16], whose dimensions can be denoted as $D_{s}$ and $D_{u}$,
respectively, then the expression

$D_{c}=D_{s}+D_{u}-N=2D_{s}-5$ (7)

is satisfied, where $D_{s}=D_{u}$ due to the symplectic nature of the
Hamiltonian under study (1), and $N=5$ is the dimension of the phase space of
the conservative system. On the other hand, if $d$ is the dimension of the
exit basin boundaries, and taking into account that such boundaries can be
described as the intersection between the stable manifold of the chaotic
saddle and the proposed initial condition plane, whose dimension is $D=2$,
then we have

$d=D+D_{s}-N=2+D_{s}-5.$ (8)

Finally, we relate the dimensions of the chaotic saddle and the basin boundary
by operating Eqs. (7) and (8), yielding

$D_{c}=2d+1.$ (9)

The latter relation enables the dimension $D_{c}$ to be computed directly if
$d$ is known, which is obtained by the uncertainty exponent algorithm as
follows [30]. We compute the exit given an arbitrarily initial condition
$(x_{0},y_{0})$ and also given four additional trajectories starting from
slightly perturbed initial conditions, such as $(x_{0}+\varepsilon,y_{0})$,
$(x_{0}-\varepsilon,y_{0})$, $(x_{0},y_{0}+\varepsilon)$ and
$(x_{0},y_{0}-\varepsilon)$, where $\varepsilon\ll 1$. If all five exits
coincide, we denote the initial condition as certain. Otherwise, we say the
initial condition is uncertain. Repeating this process for many initial
conditions and several values of the perturbation, we obtain that the fraction
of uncertain initial conditions behaves as
$f(\epsilon)\sim\varepsilon^{D-d}=\varepsilon^{2-d}$. Taking decimal
logarithms, we obtain $\log_{10}f(\epsilon)=(2-d)\log_{10}\varepsilon+k$,
where $k$ is a constant. This formula allows us to compute $d$ by a linear fit
of $\log_{10}f(\varepsilon)$ versus $\log_{10}\varepsilon$.

Two extreme situations can be distinguished here. The unpredictability is
maximal when $d=D=2$, since the basin boundary densely occupies the entire
space of initial conditions, and even under minimal system perturbations, the
fraction of uncertain conditions does not decrease, i.e., $f(\varepsilon)\sim
1$. In this situation, $D_{c}=N=5$ holds, implying that the chaotic saddle
densely occupies the entire phase space. On the other hand, the
unpredictability is minimal when the basin boundary is completely smooth,
$d=1$, and consequently $D_{c}=3$. These two extreme cases have not been
observed in any of the simulations carried out, so it is expected that the
values of the fractal dimensions are $1<d<2$, and hence $3<D_{c}<5$. This
result is in agreement with those displayed in Fig. 4.

For the sake of clarity, we only use the basin region
$x_{0},y_{0}\in[-0.1,0.1]$ to speed up the uncertainty exponent algorithm.
Specifically, for each pair of values $E_{0}$ and $\mu$, we simulate as many
initial conditions as necessary to reach one hundred uncertain conditions and,
in addition, repeat the process for fifteen values of the perturbation
$\varepsilon\in[10^{-4},10^{-1}]$.

## References

## References

  * [1] J. M. Seoane and M. A. F. Sanjuán, New developments in classical chaotic scattering, Rep. Prog. Phys. 76, 016001 (2013). doi:10.1088/0034-4885/76/1/016001.
  * [2] A. Daitche and T. Tél, Memory effects in chaotic advection of inertial particles, New J. Phys. 16, 073008 (2014). doi:10.1088/1367-2630/16/7/073008.
  * [3] J. Q. Toledo-Marín and G. G. Naumis, Escape time, relaxation, and sticky states of a softened Hénon-Heiles model: Low-frequency vibrational mode effects and glass relaxation, Phys. Rev. E 97, 042106 (2018). doi:10.1103/PhysRevE.97.042106
  * [4] I. Scheuring, T. Czárán, P. Szabó, G. Károlyi, and Z. Toroczkai, Spatial models of prebiotic evolution: soup before pizza?, Orig. Life Evol. Biosph. 33, 319–355 (2003). doi:10.1023/a:1025742505324.
  * [5] A. B. Schelin, G. Károlyi, A. P. S. de Moura, N. A. Booth, and C. Grebogi, Fractal structures in stenoses and aneurysms in blood vessels, Phil. Trans. R. Soc. A. 368, 5605–5617 (2010). doi:10.1098/rsta.2010.0268.
  * [6] T. Tél, The joy of transient chaos, Chaos 25, 097619 (2015). doi:10.1063/1.4917287.
  * [7] E. E. Zotos, Crash test for the Copenhagen problem with oblateness, Celest. Mech. Dyn. Astron. 122, 75–99 (2015). doi:10.1007/s10569-015-9611-x.
  * [8] J. F. Navarro, Dependence of the escape from an axially symmetric galaxy on the energy, Sci. Rep. 11, 8427 (2021). doi:10.1038/s41598-021-87670-5.
  * [9] A. R. Nieto, J. M. Seoane, J. E. Alvarellos, and M. A. F. Sanjuán, Resonant behavior and unpredictability in forced chaotic scattering, Phys. Rev. E 98, 062206 (2018). doi:10.1103/PhysRevE.98.062206.
  * [10] A. R. Nieto, J. M. Seoane, and M. A. F. Sanjuán, Trapping enhanced by noise in nonhyperbolic and hyperbolic chaotic scattering, Commun. Nonlinear Sci. Numer. Simulat. 102, 105905 (2021). doi:10.1016/j.cnsns.2021.105905.
  * [11] A. E. Motter and Y.-C. Lai, Dissipative chaotic scattering, Phys. Rev. E 65, 015205(R) (2001). doi:10.1103/PhysRevE.65.015205.
  * [12] J. M. Seoane, J. Aguirre, and M. A.F. Sanjuán, Basin topology in dissipative chaotic scattering, Chaos 16, 023101 (2006). doi:10.1063/1.2173342.
  * [13] J. M. Seoane, M. A. F. Sanjuán, and Y.-C. Lai, Fractal dimension in dissipative chaotic scattering, Phys. Rev. E 76, 016208 (2007). doi:10.1103/PhysRevE.76.016208.
  * [14] J. Aguirre, J. C. Vallejo, and M. A. F. Sanjuán, Wada basins and chaotic invariant sets in the Hénon-Heiles system, Phys. Rev. E 64, 066208 (2001). doi:10.1103/PhysRevE.64.066208.
  * [15] K. Yoneyama, Theory of continuous sets of points, Tokohu Math. J. 11, 43–158 (1917).
  * [16] J. Aguirre, R. L. Viana, and M. A. F. Sanjuán, Fractal structures in nonlinear dynamics, Rev. Mod. Phys. 81, 333 (2009). doi:10.1103/RevModPhys.81.333.
  * [17] A. Wagemakers, A. Daza, and M. A.F. Sanjuán, How to detect Wada basins, Discrete Continuous Dyn. Syst. Ser. B 26(1), 717–739 (2021). doi:10.3934/dcdsb.2020330
  * [18] Z. Kovács and L. Wiesenfeld, Topological aspects of chaotic scattering in higher dimensions, Phys. Rev. E 63, 056207 (2001). doi:10.1103/PhysRevE.63.056207.
  * [19] Y.-D. Lin, A. M. Barr, L. E. Reichl, and C. Jung, Fractal scattering dynamics of the three-dimensional HOCl molecule, Phys. Rev. E 87, 012917 (2013). doi:10.1103/PhysRevE.87.012917.
  * [20] G. Drótos, F. G. Montoya, C. Jung, and T. Tél, Asymptotic observability of low-dimensional powder chaos in a three-degrees-of-freedom scattering system, Phys. Rev. E 90, 022906 (2014). doi:10.1103/PhysRevE.90.022906.
  * [21] S. Naik and S. Wiggins, Finding normally hyperbolic invariant manifolds in two and three degrees of freedom with Hénon-Heiles-type potential, Phys. Rev. E 100, 022204 (2019). doi:10.1103/PhysRevE.100.022204.
  * [22] Y.-C. Lai, A. P. S. de Moura, and C. Grebogi, Topology of high-dimensional chaotic scattering, Phys. Rev. E 62, 6421 (2000). doi:10.1103/PhysRevE.62.6421.
  * [23] A. P. S. de Moura and C. Grebogi, Countable and uncountable boundaries in chaotic scattering, Phys. Rev. E 66, 046214 (2002). doi:10.1103/PhysRevE.66.046214.
  * [24] D. ter Haar, The Vibrational Levels of an Anharmonic Oscillator, Phys. Rev. 70, 222 (1946). doi:10.1103/PhysRev.70.222.
  * [25] D. S. Fernández, A. G. López, J. M. Seoane, and M. A.F. Sanjuán, Ergodic decay laws in Newtonian and relativistic chaotic scattering, Commun. Nonlinear Sci. Numer. Simulat. 103, 105987 (2021). doi:10.1016/j.cnsns.2021.105987.
  * [26] H. E. Kandrup and S. J. Novotny, Phase Mixing in Unperturbed and Perturbed Hamiltonian Systems, Celest. Mech. Dyn. Astron. 88, 1 (2004). doi:10.1023/B:CELE.0000009380.17257.98.
  * [27] R. L. Burden, J. D. Faires, and A. M. Burden, Numerical Analysis, Tenth Edition, Cengage Learning (2015). doi:10.13140/2.1.4830.2406.
  * [28] E. Ott, _Chaos in Dynamical Systems_ , Second Edition, Cambridge University Press (2002). doi:10.1017/CBO9780511803260.
  * [29] A. R. Nieto, E. E. Zotos, J. M. Seoane, and M. A. F. Sanjuán, Measuring the transition between nonhyperbolic and hyperbolic regimes in open Hamiltonian systems, Nonlinear Dyn. 99, 3029–3039 (2020). doi:10.1007/s11071-019-05433-4.
  * [30] S. W. McDonald, C. Grebogi, E. Ott, and J. A. Yorke, Fractal basin boundaries, Physica D 17, 125 (1985). doi:10.1016/0167-2789(85)90001-6.
  * [31] Y. Zhang and G. Luo, Wada bifurcations and partially Wada basin boundaries in a two-dimensional cubic map, Phys. Lett. A 377, 1274–1281 (2013). doi:10.1016/j.physleta.2013.03.027.
  * [32] A. Daza, J. O. Shipley, S. R. Dolan, and M. A. F. Sanjuán, Wada structures in a binary black hole system, Phys. Rev. D 98, 084050 (2018). doi:10.1103/PhysRevD.98.084050.
  * [33] X. Chen, T. Nishikawa, and A. E. Motter, Slim Fractals: The Geometry of Doubly Transient Chaos, Phys. Rev. X 7, 021040 (2017). doi:10.1103/PhysRevX.7.021040.
  * [34] A. P. S. de Moura and C. Grebogi, Chemical and biological activity in three-dimensional flows, Phys. Rev. E 70, 026218 (2004). doi:10.1103/PhysRevE.70.026218.
  * [35] R. Klages, S. S. G. Gallegos, J. Solanpaa, M. Sarvilahti, and E. Rasanen, Normal and Anomalous Diffusion in Soft Lorentz Gases, Phys. Rev. Lett. 122, 064102 (2019). doi:10.1103/PhysRevLett.122.064102.
  * [36] A. Daza, A. Wagemakers, and M. A. F. Sanjuán, Wada property in systems with delay, Commun. Nonlinear Sci. Numer. Simulat. 43, 220–226 (2017). doi:10.1016/j.cnsns.2016.07.008.


Hitting Time Distribution for finite states Markov Chain111The project is
partially supported by the National Natural Science Foundation of China (Grant
No. 11131003).

Wenming Hong222School of Mathematical Sciences & Laboratory of Mathematics and
Complex Systems, Beijing Normal University, Beijing 100875, P.R. China. Email:
<EMAIL_ADDRESS>Ke Zhou333 School of Mathematical Sciences & Laboratory of
Mathematics and Complex Systems, Beijing Normal University, Beijing 100875,
P.R. China<EMAIL_ADDRESS>
Abstract

Consider a Markov chain with finite state $\\{0,1,\cdots,d\\}$ and absorbing
at state $d$. We give out the generation functions (or Laplace transforms) for
the absorbing time when the chain starts from any state $i$. The results
generalize the well-known theorems for the birth-death (Karlin and McGregor
[7], 1959) and the skip-free (Brown and Shao [1], 1987) Markov chain starts
from state $0$. Our proof is directly and simple.

Keywords: Markov chain, absorbing time, generation functions, Laplace
transforms, eigenvalues.

Mathematics Subject Classification (2010): 60E10, 60J10, 60J27, 60J35.

##  1  Introduction

For the birth-death and the skip-free (upward jumps may be only of unit size,
and there is no restriction on downward jumps) Markov chain with finite state
$\\{0,1,\cdots,d\\}$ and absorbing at state $d$, an well-known interesting
property for the absorbing time is that it is distributed as a summation of
$d$ independent geometric (or exponential) random variables.

There are many authors give out different proofs to the results. For the birth
and death chain, the well-known results can be traced back to Karlin and
McGregor ([7], 1959) Keilson ([8], 1971; [9]). Kent and Longford ([10], 1983)
proved the result for the discrete time version (nearest random walk) although
they have not specified the result as usual form (section 2, [10]). Fill ([4],
2009) gave the first stochastic proof to both nearest random walk and birth
and death chain cases via duality which was established in [2]. Diaconis and
Miclo ([3], 2009) presented another probabilistic proof for birth and death
chain. Very recently, Gong, Mao and Zhang ([6], 2012) gave a similar result in
the case that the state space is $\mathbb{Z^{+}}$. For the skip-free chain,
Brown and Shao ([1], 1987) first proved the result in continuous time
situation; Fill ([5], 2009) gave a stochastic proof to both discrete and
continuous time cases also by using the duality, and considered the general
finite-state Markov chain situation when the chain starts from the state $0$.

The purpose of this paper is to consider a general Markov chain with finite-
state $\\{0,1,\cdots,d\\}$ and absorbing at state $d$. We give out the
generation functions (or Laplace transforms) for the absorbing time when the
chain starts from any state $i$ (not just from state $0$ only). In
particulary, the results generalize the well-known theorems for the birth-
death (Karlin and McGregor [7], 1959) and the skip-free (Brown and Shao [1],
1987) Markov chain.

Our proof is to calculate directly the generation functions (or Laplace
transforms) of the absorbing times by iterating method, which have been used
for the skip-free Markov chain by Zhou ([11]). After revised the method in
[11], we found the proof is very simple and enable us to deal with the general
finite-state Markov chain starts from any state $i$. The key idea is to
consider directly the absorbing time $\tau_{i,d}$ starting from any state $i$.

##  2  Discrete time

Define the transition probability matrix $P$ as

$P=\left(\begin{array}[]{ccccccc}r_{0}&p_{0,1}&p_{0,2}&&\cdots&p_{0,d-1}&p_{0,d}\\\
q_{1,0}&r_{1}&p_{1,2}&&\cdots&p_{1,d-1}&p_{1,d}\\\
\vdots&\ddots&\ddots&&\ddots&\ddots&\vdots\\\
q_{d-1,0}&q_{d-1,1}&q_{d-1,2}&&\cdots&r_{d-1}&p_{d-1,d}\\\
0&0&0&&\cdots&0&1\\\ \end{array}\right)_{(d+1)\times(d+1)},$ (2.1)

and for $1\leq j\leq d+1$, we denote $A_{j}(s)$ as the $d\times d$ sub-matrix
by deleting the $(d+1)^{th}$ line and the $j^{th}$ row of the matrix
$I_{d+1}-sP$.

Let $\tau_{i,d}$ be the absorbing time of state $d$ starting from $i$ and
$f_{i}(s)$ be the generation function of $\tau_{i,d}$,

$f_{i}(s)=\mathbb{E}s^{\tau_{i,d}}~{}~{}~{}\text{for}~{}0\leq i\leq d-1.$
(2.2)

######  Theorem 2.1.

For $0\leq i\leq d-1$,

$f_{i}(s)=(-1)^{d-i}\frac{\det A_{i+1}(s)}{\det A_{d+1}(s)}.$ (2.3)

###### Proof.

By decomposing the first step ,the generation function of $\tau_{i,d}$
satisfy, for $0\leq i\leq d-1$,

$\begin{split}f_{i}(s)&=r_{i}sf_{i}(s)+p_{i,i+1}sf_{i+1}(s)+p_{i,i+2}sf_{i+2}(s)+\cdots
p_{i,d-1}sf_{d-1}(s)+p_{i,d}s\\\
&~{}~{}~{}~{}~{}q_{i,i-1}sf_{i-1}(s)+q_{i,i-2}sf_{i-2}(s)+\cdots+q_{i,0}sf_{0}(s).\end{split}$
(2.4)

These system of $d$ equations are linear about
$f_{0}(s),~{}f_{1}(s)\cdots,~{}f_{d-1}(s)$. Use Cramer’s Rule, we can solve
from these system of $d$ equations and get (2.3). $\Box$

Immediately, we obtain the results for the skip-free discrete time Markov
chain (Fill [5], 2009).

######  Corollary 2.1.

If we assume for $j-i>1$, $p_{i,j}=0$. We have

$f_{0}(s)=\prod_{i=0}^{d-1}\left[\frac{(1-\lambda_{i})s}{1-\lambda_{i}s}\right],$
(2.5)

where $\lambda_{0},\cdots,\lambda_{d-1}$ are the $d$ non-unit eigenvalues of
$P$.

In particular, if all of the eigenvalues are real and nonnegative, then the
hitting time is distributed as the sum of $d$ independent geometric random
variables with parameters $1-\lambda_{i}$.

###### Proof.

Note that, $1$ is evidently an eigenvalue of $P$. So on the one hand
${\det(I_{d+1}-sP)}=(1-s)\prod_{i=0}^{d-1}(1-\lambda_{i}s)$ (where
$\lambda_{0},\cdots,\lambda_{d-1}$ are the $d$ non-unit eigenvalues of $P$);
on the other hand from (2.1) we have ${\det(I_{d+1}-sP)}=(1-s)\times\det
A_{d+1}(s)$; as a consequence we get

${\det A_{d+1}(s)}=\prod_{i=0}^{d-1}(1-\lambda_{i}s).$ (2.6)

From (2.1) and the definition of $A_{j}$, it is easy to get

$\det A_{1}(s)=(-1)^{d}p_{0,1}p_{1,2}\cdots p_{d-1,d}s^{d}.$ (2.7)

By (2.3) and (2.6) we have

$\displaystyle\det A_{1}(1)=(-1)^{d}f_{0}(1)\cdot{\det
A_{d+1}(1)}=(-1)^{d}f_{0}(1)\cdot\prod_{i=0}^{d-1}(1-\lambda_{i}).$

On the other hand, from (2.7) we get $\det
A_{1}(1)=(-1)^{d}p_{0,1}p_{1,2}\cdots p_{d-1,d}$ and recall that $f_{0}(1)=1$
by (2.2), we obtain that

$p_{0,1}p_{1,2}\cdots p_{d-1,d}=\prod_{i=0}^{d-1}(1-\lambda_{i}).$ (2.8)

Then by (2.7) and (2.8)

$\text{det}A_{1}(s)=(-1)^{d}\prod_{i=0}^{d-1}(1-\lambda_{i})s^{d},$ (2.9)

and (2.5) is immediate from (2.6) and (2.9). $\Box$

##  3  Continuous time

Define the generator $Q$ as

$Q=\left(\begin{array}[]{ccccccc}-\gamma_{0}&\alpha_{0,1}&\alpha_{0,2}&&\cdots&\alpha_{0,d-1}&\alpha_{0,d}\\\
\beta_{1,0}&-\gamma_{1}&\alpha_{1,2}&&\cdots&\alpha_{1,d-1}&\alpha_{1,d}\\\
\vdots&\ddots&\ddots&&\ddots&\ddots&\vdots\\\
\beta_{d-1,0}&\beta_{d-1,1}&\beta_{d-1,2}&&\cdots&-\gamma_{d-1}&\alpha_{d-1,d}\\\
0&0&0&&\cdots&0&0\\\ \end{array}\right)_{(d+1)\times(d+1)},$

and for $1\leq j\leq d+1$, we denote $\widetilde{A}_{j}(s)$ as the $d\times d$
sub-matrix by deleting the $(d+1)^{th}$ line and the $j^{th}$ row of the
matrix $I_{d+1}-sP$. Let $\tau_{i,d}$ be the absorbing time of state $d$
starting from $i$ and $\widetilde{f}_{i}(s)$ be the Laplace transform of
$\tau_{i,d}$.

$\widetilde{f}_{i}(s)=\mathbb{E}e^{-s\tau_{i,d}}.$

It is well known that the chain on the finite state has an simple structure.
The process starts at $i$, it stay there with an
$\mbox{Exponential}~{}(\gamma_{i})$ time, then jumps to $i+j$ with probability
$\frac{\alpha_{i,i+j}}{\gamma_{i}}$, to $i-k$ with probability
$\frac{\beta_{i,i-k}}{\gamma_{i}}$ .

######  Theorem 3.1.

$\widetilde{f}_{i}(s)=(-1)^{d-i}\frac{\det\widetilde{A}_{i+1}}{\det\widetilde{A}_{d+1}},~{}~{}~{}\text{for}~{}0\leq
i\leq d-1,$ (3.1)

###### Proof.

By decomposing the trajectory at the first jump, $0\leq i\leq d-1$,

$\begin{split}\widetilde{f}_{i}(s)&=\frac{\gamma_{i}}{\gamma_{i}+s}\frac{\alpha_{i,i+1}}{\gamma_{i}}\widetilde{f}_{i+1}(s)+\frac{\gamma_{i}}{\gamma_{i}+s}\frac{\alpha_{i,i+2}}{\gamma_{i}}\widetilde{f}_{i+1}(s)+\cdots+\frac{\gamma_{i}}{\gamma_{i}+s}\frac{\alpha_{i,d-1}}{\gamma_{i}}\widetilde{f}_{d-1}(s)+\frac{\gamma_{i}}{\gamma_{i}+s}\frac{\alpha_{i,d}}{\gamma_{i}}+\\\
&~{}~{}~{}~{}~{}+\frac{\gamma_{i}}{\gamma_{i}+s}\frac{\beta_{i,i-1}}{\gamma_{i}}\widetilde{f}_{i-1}(s)+\frac{\gamma_{i}}{\gamma_{i}+s}\frac{\beta_{i,i-2}}{\gamma_{i}}\widetilde{f}_{i-2}(s)+\cdots+\frac{\gamma_{i}}{\gamma_{i}+s}\frac{\beta_{i,0}}{\gamma_{i}}\widetilde{f}_{0}(s)\\\
&=\frac{\alpha_{i,i+1}}{\gamma_{i}+s}\widetilde{f}_{i+1}(s)+\cdots+\frac{\alpha_{i,d-1}}{\gamma_{i}+s}\widetilde{f}_{d-1}(s)+\frac{\alpha_{i,d}}{\gamma_{i}+s}+\frac{\beta_{i,i-1}}{\gamma_{i}+s}\widetilde{f}_{i-1}(s)+\cdots+\frac{\beta_{i,0}}{\gamma_{i}+s}\widetilde{f}_{0}(s)\end{split}$

These system of $d$ equations are linear about
$\widetilde{f}_{0}(s),~{}\widetilde{f}_{1}(s)\cdots,~{}\widetilde{f}_{d-1}(s)$.
Use Cramer’s Rule, we can solve from these system of $d$ equations and get
(3.1). $\Box$

Immediately, we obtain the results for the skip-free continuous time Markov
chain (Brown and Shao [1], 1987).

######  Corollary 3.1.

If we assume for $j-i>1$, $\alpha_{i,j}=0$. We have

$\varphi_{d}(s)=\prod_{i=0}^{d-1}\frac{\lambda_{i}}{s-\lambda_{i}},$

where $\lambda_{i}$ are the $d$ non-zero eigenvalues of $-Q$.

In particular, if all of the eigenvalues are real and nonnegative, then the
hitting time is distributed as the sum of $d$ independent exponential random
variables with parameters $\lambda_{i}$.

###### Proof.

The proof is similar as Corollary 2.1, we can calculate that
${\det\widetilde{A}_{d+1}}=\prod_{i=0}^{d-1}(s-\lambda_{i})$,
$\det\widetilde{A}_{1}=(-1)^{d}\alpha_{0,1}\alpha_{1,2}\cdots\alpha_{d-1,d}=(-1)^{d}\prod_{i=0}^{d-1}\lambda_{i}$.
$\Box$

## References

  * [1] Brown, M. and Shao, Y. S. Identifying coefficients in the spectral representation for first passage time distributions. _Probab. Eng. Inform. Sci._ 1 (1987), 69–74.
  * [2] Diaconis, P. and Fill, J. A. Strong stationary times via a new form of duality. _Ann. Probab._ 18 (1990), 1483–1522.
  * [3] Diaconis, P. and Miclo, L. On times to quasi-stationarity for birth and death processes. _J. Theoret. Probab._ 22 (2009), 558–586.
  * [4] Fill, J. A. The passage time distribution for a birth-and-death chain: Strong stationary duality gives a first stochastic proof. _J. Theoret. Probab._ 22 (2009), 543–557.
  * [5] Fill, J. A. On hitting times and fastest strong stationary times for skip-free and more general chains. _J. Theoret. Probab._ 22 (2009), 587–600.
  * [6] Gong, Y. Mao, Y. H and Zhang, C. Hitting time distributions for denumerable birth and death processes. _J. Theoret. Probab._ 25 (2012), 950–980.
  * [7] Karlin, S. and McGregor, J. Coincidence properties of birth and death processes. _Pacific J. Math._ 9 (1959), 1109–1140.
  * [8] Keilson, J. Log-concavity and log-convexity in passage time densities for of diffusion and birth-death processes. _J. Appl. Probab._ 8 (1971), 391–398.
  * [9] Keilson, J.  _Markov Chain Models—Rarity and Exponentiality._ Springer, New York, 1979.
  * [10] Kent, J. T. and Longford, N. T An eigenvalue decomposition for first hitting times in random walk. _Z. Wahrscheinlichkeitstheorie verw. Gebiete_ 63 (1983), 71–84.
  * [11] Zhou, K (2013), Hitting Time Distribution for Skip-Free Markov Chains: A Simple Proof. _submitted._


# Integrating LSTM and BERT for Long-Sequence Data Analysis in Intelligent
Tutoring Systems ††thanks: Citation: Authors. Title. Pages…. DOI:000000/11111.

Zhaoxing Li, Sebastian Stein
School of Electronics and Computer Science
University of Southampton
Southampton
{zhaoxing.li<EMAIL_ADDRESS>
&Jujie Yang, Jindi Wang
Department of Computer Science
Durham University
Durham
{jujie.yang<EMAIL_ADDRESS>
Lei Shi
School of Computing
Newcastle University
Newcastle upon Tyne
<EMAIL_ADDRESS>

###### Abstract

The field of Knowledge Tracing aims to understand how students learn and
master knowledge over time by analyzing their historical behaviour data. To
achieve this goal, many researchers have proposed Knowledge Tracing models
that use data from Intelligent Tutoring Systems to predict students’
subsequent actions. However, with the development of Intelligent Tutoring
Systems, large-scale datasets containing long-sequence data began to emerge.
Recent deep learning based Knowledge Tracing models face obstacles such as low
efficiency, low accuracy, and low interpretability when dealing with large-
scale datasets containing long-sequence data. To address these issues and
promote the sustainable development of Intelligent Tutoring Systems, we
propose a LSTM BERT-based Knowledge Tracing model for long sequence data
processing, namely LBKT, which uses a BERT-based architecture with a Rasch
model-based embeddings block to deal with different difficulty levels
information and an LSTM block to process the sequential characteristic in
students’ actions. LBKT achieves the best performance on most benchmark
datasets on the metrics of ACC and AUC. Additionally, an ablation study is
conducted to analyse the impact of each component of LBKT’s overall
performance. Moreover, we used t-SNE as the visualisation tool to demonstrate
the model’s embedding strategy. The results indicate that LBKT is faster, more
interpretable, and has a lower memory cost than the traditional deep learning
based Knowledge Tracing methods.

_K_ eywords Knowledge Tracing $\cdot$ BERT $\cdot$ Student Modelling $\cdot$
Long-Sequence Data Processing $\cdot$ Technology Enhanced Learning (TEL)

##  1 Introduction

Technology Enhanced Learning (TEL) has become increasingly important in
providing high-quality education to build a more sustainable world. The recent
COVID-19 pandemic has significantly impacted traditional classroom education
and sparked online learning, enabling teaching and learning remotely.
Meanwhile, the development of online learning systems has made it possible to
use Intelligent Tutoring Systems (ITS) to store and analyse a sizable amount
of student behaviour data to improve intelligent educational services. As one
of the widely applied TEL technologies, Knowledge Tracing (KT) has drawn a lot
of attention. KT is the field of modelling students’ learning trajectories and
predicting their sequential actions based on historical interaction data
between students and ITS [1].

With the development of ITS, large-scale datasets such as _EdNet_ [2] and
_Junyi Academy_ [3] began to emerge. In these datasets, long-sequence student
interaction data were gathered as an increasing number of students used the
ITS for an extended period. The long- and short-sequence data in these
datasets are unbalanced, which satisfies the long-tail distribution [4]. For
instance, within the EdNet dataset, a substantial amount of student action
sequences are included, ranging from the shortest sequence that may comprise
just a single action to the longest sequence that encompasses 40,157 actions.
Notably, the average action sequence length of the EdNet dataset is 121.5,
indicating a moderate length of data sequences overall. However, it is
important to note that the distribution of sequence lengths is highly skewed,
and this unbalanced distribution has an impact on the overall performance of
the KT models [5]. Although the quantity of short-sequence data is larger than
the long-sequence data, the latter is of more weight than the former in
prediction tasks [4].

In general, KT models could be divided into three categories: probabilistic KT
models, logistic KT models, and deep learning based KT methods (DKT) [6].
Traditional probabilistic KT models and logistic KT models are forced to
confront difficulties such as decreased processing efficiency and increased
memory usage as growing amounts of longer sequence data are released. Deep
learning based KT models are known to suffer from inefficiencies when
processing long-sequence action data problems, including issues related to the
accuracy, speed, and memory usage [5, 4]. Therefore, allowing the processing
of very long sequence data is key to achieving high performance for next-
generation KT models [5]. Moreover, due to the black-box nature of traditional
deep learning methods, the current deep learning based KT models also struggle
with the lack of interpretability [7].

To address the above issues, in this paper, we propose LBKT, a novel LSTM BERT
Knowledge Tracing model, for processing long sequence data. The model combines
the strength of the Bidirectional Encoder Representations from Transformers
(BERT) model in capturing the relations of complex data [8] with the strength
of the LSTM model in handling long sequential data to improve its performance
on large-scale datasets containing long-sequence data (here, the long-sequence
data indicates a length longer than 400 interactions). Moreover, we utilise a
Rasch model-based embedding method to process the difficulty level information
in the historical behaviour data of students. The Rasch model is a classic yet
powerful model in psychometrics [9], which could be utilised to construct raw
questions and knowledge embeddings for KT tasks [7]. Rasch model based
embedding could improve the model’s performance and interpretability. The
experimental results show that our proposed LBKT outperforms the baseline
models in five datasets on metrics ACC and AUC. Moreover, it is faster at
processing long-sequence data at two long-sequence datasets we extract from
the two large-scale datasets. Furthermore, we use t-SNE as the visualisation
tool to demonstrate the interpretability of the embedding strategy.

The main contributions of our paper lie in the following two aspects:

  1. 1.

We propose LBKT 111Source code and datasets are available at
https://github.com/******/LBKT, a novel LSTM BERT Knowledge Tracing model for
long sequence data processing. The LBKT leverages the power of BERT, Rasch-
based embedding strategies, and LSTM.

  2. 2.

The experimental results show that LBKT outperforms the baseline models on
five ITS datasets on the metric of AUC(assist12, assist17, algebra06, EdNet,
and Junyi Academy). Another comparative experiments show the effectiveness of
LBKT when processing long-sequence datasets. LBKT model exhibits better
interpretability than traditional deep learning based KT models and has
advantages in training efficiency.

##  2 Related Work

###  2.1 Knowledge Tracing

Knowledge Tracing (KT) is used in Intelligent Tutoring Systems (ITS) to model
and predict a student’s mastery level of a specific skill or concept over time
[10]. It is based on the assumption that a student’s knowledge state is a
hidden variable that can be inferred from their observable behaviour, such as
their responses to questions or tasks related to the skill or concept being
measured [11]. Its goal is to provide personalised feedback and support to
students by tracking their progress and adapting instruction to meet their
individual needs. This can help to improve student learning outcomes and
enhance educational effectiveness. Broadly, there are three categories of KT
methods: probabilistic KT models, logistic KT models, and deep learning-based
KT models [6].

Probabilistic KT models assume that the student’s learning process follows a
Markov Process, where students’ knowledge mastery could be measured by their
observed learning performance [11]. Bayesian KT, or BKT, is the earliest and
most classic probabilistic model, which was inspired by cognitive mastery
learning [12]. BKT models generally use a probabilistic graphical model, such
as Hidden Markov Model (HMM)[11] and Bayesian Belief Network [13], to track
students’ changing learning states. The major shortcoming of BKT is that it
assumes a simplistic two-state student modelling framework, where a student’s
knowledge is either learned or unlearned, and there is no concept of
forgetting or decay in the model. However, in reality, a student’s knowledge
could be complex and multi-faceted and could change over time due to various
factors such as decay and interference. Therefore, BKT may not be able to
capture the nuances of student learning and may not provide an accurate
representation of their knowledge state over time. For example, BKT assumes
that each question only required one skill and that the various skills were
irrelevant to each other [11, 14]. Therefore, in general, BKT models cannot
process complicated problems, including the multiple skills and the complex
relationship among the concepts, questions, and skills. To address this
limitation, Käser et al. proposed Dynamic BKT, or DBKT, based on Dynamic
Bayesian Network (DBN), to model the prerequisite hierarchies and dependencies
of multiple skills [15]. However, both BKT and DBKT still struggle with
processing multiple topics or skills, failing to account for contextual
factors that may impact student learning.

The logistic KT models are built on the principle of logistic regression,
which is a statistical method used to model the probability of a binary
outcome based on one or more predictor variables [6]. In the context of
educational data, the predictor variables could include a student’s prior
performance on a set of related skills or concepts, their response time, and
their correctness or incorrectness in answering assessment questions. The
output of the logistic regression KT model is a probability estimate of a
student’s mastery level on a particular skill or concept, which can be used to
inform personalized learning interventions and improve student outcomes. There
are three logistic models. The Learning Factor Analysis model (LFA)
incorporates the initial knowledge state, easiness of knowledge components
(KCs), and learning rate of KCs to estimate the student’s initial knowledge
state, the easiness of different KCs, and the learning rate of KCs [16]. The
Performance Factor Analysis (PFA) model is an extension of the LFA model and
takes into account the student’s performance. PFA considers parameters for
previous failures (f) and successes (s) for the KC, in addition to the
easiness of KCs [17]. The Knowledge Tracing Machines (KTM) model uses
factorization machines (FMs) to extend logistic models to higher dimensions
[15].

Inspired by the recent success of deep learning (DL)[18], researchers have
applied deep learning technologies into the KT field to develop DL-based
Knowledge Tracing [19]. DL-based KT typically models a knowledge tracing task
as a sequence prediction problem. With the self-attention architectures
applied in the deep learning field, KT models based on the self-attention
mechanism began to emerge. For example, SAKT [20] and SAINT+ [21] apply the
self-attention mechanism to KT models and achieve higher performance than the
traditional DL-based methods. With the development of the self-attention
mechanism, Transformer based knowledge tracing models also have been proposed.
Ghosh[7] proposes context-aware attentive knowledge tracing (AKT), which
introduces a novel monotonic attention mechanism that accounts for the
temporal nature of the learning process and the decay of students’ knowledge.
Nakagawa et al. proposed the Graph-based Knowledge Tracing (GKT) model, which
incorporates the potential graph structure of KCs into a graph [22]. There
were also KT methods based on BERT that had been proposed. MonacoBERT [23] is
a BERT-based KT model that incorporates the monotonic convolutional multi-head
attention and classical test-theory-based (CTT-based) embedding strategy to
improve performance. BEKT [24] is a Bidirectional Encoder representation from
the Transformers-based model that predicts student knowledge state by
combining historical learning performance.

###  2.2 Transformer-based Model and Application

Transformer is a prominent neural network model proposed by Vaswani et al.,
which utilises the self-attention mechanism to extract inherent features.
Transformer-based models have achieved significant success in the Deep
Learning field, especially in Nature Language Processing (NLP) and image
generation tasks [25, 26].

The evolution of Transformer-based models, such as BERT [8] and GPT[27], has
achieved outstanding performance in the above tasks. BERT, first proposed by
Devlin et al., is a successful application of Transformer [8]. BERT utilises
the self-attention mechanism and the masked language model (MLM) to train the
Transformer bidirectionally in the NLP fields [8]. BERT is renowned for its
exceptional ability to process and comprehend natural language text
efficiently. It has consistently outperformed other deep learning models in a
broad range of tasks, extending beyond the field of NLP. BERT’s success can be
attributed to several key features, including its bidirectional context, which
allows it to capture the dependencies between both preceding and succeeding
tokens in a sequence. Additionally, BERT’s large pre-training corpus enables
it to learn a robust language representation that can be fine-tuned for
downstream tasks with relatively small amounts of labelled data. BERT’s
transformer architecture, which uses self-attention mechanisms to capture
global dependencies between tokens in a sequence, is also a significant factor
contributing to its performance. The self-attention mechanism allows BERT to
weigh the importance of different tokens in a sequence dynamically, which
improves its ability to capture complex patterns and relationships in the data
[8]. BERT is also known for its ability to generate high-quality embeddings,
which are crucial for many natural language processing tasks [8]. There have
been a lot of BERT variants applied in other deep learning fields,
demonstrating their outstanding performances. For example, ConvBERT [28]
applies the original BERT architecture in the image processing field; BERT4Rec
uses BERT model to improve recommendation systems [29]; LakhNES uses BERT
model to enhance Music Generation [30]. However, in the Knowledge Tracing
field, although some BERT-based models, such as BEKT [31, 32, 33, 34, 35, 36,
37, 38, 39] and BiDKT [40], are proposed to improve performance, they are
unable to outperform state-of-the-art KT methods in large-scale datasets
containing long-sequence data.

##  3 Methodology

###  3.1 Problem Statement

The key to knowledge tracing is to predict the correctness of a student’s next
answer in a sequence. Let $x_{1},\ldots,x_{t}$ denote the student’s actions,
and let the $t$-th action be represented as $x_{t}=(q_{t},a_{t})$, where
$q_{t}$ is the question presented to the student and $a_{t}$ is the student’s
response. The goal is to estimate the correctness
$P(a_{t}=1|x_{1},\ldots,x_{t-1})$, that is, the correctness of student’s
response to the current question, given their previous actions in the
sequence.

###  3.2 Proposed Model Architecture

We propose a novel model, LBKT, for the task of knowledge tracing on large-
scale datasets containing long-sequence data. While previous BERT-based KT
models have shown remarkable success in capturing the relations of complex
data, they also have inefficiencies when dealing with long sequence student
action data [24]. On the other hand, LSTM models have been proven to excel in
handling long sequential data. In response to these challenges, we propose a
novel KT model that combines the strengths of both the BERT and LSTM models to
improve performance on large-scale datasets containing long-sequence data
(where long-sequence data indicates a length longer than 400 interactions).
The Rasch embedding (also known as the 1PL IRT model) is a method to represent
questions and concepts in a mathematical space [9]. The embeddings are created
using a vector that summarizes the variation in questions covering a concept
and a scalar difficulty parameter that controls how far a question deviates
from the concept it covers. The embeddings are used as raw embeddings for
questions and responses, which is a way to track a learner’s knowledge state.
By leveraging the strengths of a BERT-based model, Rasch model-based
embeddings, and long short-term memory (LSTM) unit, our proposed model
architecture has the potential to effectively process and understand
relationships among different features in long-sequence data, as illustrated
in Fig. 1.

Figure 1: The architecture of LBKT. LBKT consists of three components: 1) the
Rasch model-based embeddings (on the left), 2) the BERT-based architecture (in
the middle), and 3) the LSTM block (on the right).

The first component of LBKT is the Rasch model-based embeddings proposed by
Ghosh [7]. The Rasch model-based embeddings consist of difficulty level
embeddings $E_{\text{d}}$ and question embeddings $E_{\text{q}}$. These
embeddings are multiplied and added to the BERT token embeddings and the $sin$
and $cos$ positional embeddings to build the final embeddings, as shown in the
following equation:

$E=E_{\text{Rasch}}+E_{\text{BERT Token}}+E_{\text{Position}}$ (1)

where the Rasch model-based embeddings $E_{\text{Rasch}}$ are defined as:

$E_{\text{Rasch}}=E_{\text{d}}+E_{\text{d}}\times E_{\text{q}}$ (2)

The segment embeddings, which are typically used to represent information
about the segment in the BERT model, are replaced by the Rasch embeddings
mentioned above in our model’s architecture. Rasch model-based embeddings are
able to more accurately estimate students’ knowledge states, as explained
earlier, making them a key contributor to the effectiveness of LBKT for
knowledge tracing tasks.

The second component of LBKT is a BERT-based block, which consists of 12
Transformer blocks. Each includes a multi-head attention mechanism, a
feedforward network (FFN), and sublayer connections. The multi-head attention
mechanism uses the “Scaled Dot Product Attention” method as implemented in
BERT, along with queries $Q$, keys $K$, values $V$, and an attention mask for
padded tokens. The FFN has a feedforward hidden layer with a size of four
times that of the model’s hidden layer and uses the GELU activation function
rather than RELU.

The sublayer connections in the Transformer block include a residual
connection followed by layer normalization. The formulas for the attention
mechanism and the FFN are as follows:

$\operatorname{Attention}(Q,K,V)=\operatorname{softmax}\left(\frac{QK^{T}}{\sqrt{d_{k}}}\right)V$
(3)

$\text{FFN}(x)=\text{GELU}(W_{1}x+b_{1})W_{2}+b_{2}$ (4)

In the third component of LBKT, we use a neural network (NN) linear
transformation instead of the attention projection typically used in
conjunction with the LSTM unit. This is based on our observed improved
performance with the NN linear transformation in our experiments. It should be
noted that this choice is not necessarily related to the length or complexity
of the sequence but rather to the specific characteristics of the data and the
task at hand.

Overall, LBKT is a model that is tailored specifically for use in the field of
knowledge tracing. It combines the natural language processing capabilities of
the BERT model with the ability to accurately estimate knowledge states using
Rasch model-based embeddings and the ability to effectively handle long
sequences of data using the LSTM unit and the NN linear transformation. This
makes it an ideal choice for the task of knowledge tracing in large-scale
datasets containing long-sequence data with unbalanced data distribution.

###  3.3 Experiment Setting

####  3.3.1 Datasets

We used five benchmark datasets to validate the effectiveness of the LBKT
model, including
assist12222https://sites.google.com/site/assistmentsdata/home,
assist17333https://sites.google.com/site/assistmentsdata/home, algebra06
444https://pslcdatashop.web.cmu.edu/KDDCup, EdNet [2]
555https://github.com/riiid/ednet, and Junyi
Academy[3]666https://pslcdatashop.web.cmu.edu/Files?datasetId=1275. Table
LABEL:dataset shows the sizes of the above datasets. In general datasets, such
as assist 12 and assist 17, it could be challenging to identify and extract
large amounts of long-sequence data. Therefore, we validated the speed
performance of every model on two datasets with long-sequence student action
data extracted from EdNet and Junyi Academy. The mean action sequence length
of EdNet is 121.5. The mean interaction length of Junyi Academic is 104.7.
Table 2 shows the action sequence length statistics of EdNet and Junyi
Academy. Here, we define the longer action sequence as longer than 100
records. We extract 200 students’ action sequences that include interactions
longer than 100 actions from each dataset as the long-sequence dataset to
validate the performance of different KT models. Lastly, we selected different
lengths of action sequences from Ednet to test the speed performance of each
model. We selected four groups with average records lengths of 100, 200, 300,
and 400, respectively. Each of these groups included 50 students.

Table 1: Benchmark dataset data statistics Dataset | Students | Concepts | Questions | Interactions
---|---|---|---|---
assist12 | 24,429 | 264 | 51,632 | 1,968,737
assist17 | 1,708 | 411 | 3,162 | 934,638
algebra06 | 1,318 | 1,575 | 549,821 | 1,808,533
EdNet | 784,309 | 1,472 | 11,957 | 641,712
Junyi Academy | 247,606 | 13,169 | 722 | 25,925,922
Table 2: Action Length Statistics of EdNet and Junyi Academy Features | Junyi Academy | EdNet
---|---|---
Students | 247606 | 784309
Max action length | 22067 | 40157
Mean action length | 104.7 | 121.5

####  3.3.2 Baseline Models

We compared our LBKT to three state-of-the-art models, BEKT [24], AKT [7],
DKVMN [41], as well as the two top baseline models in the Riiid Answer
Correctness Prediction Competition provided by
Kaggle777https://www.kaggle.com/code/datakite/riiid-answer-correctness,
including SSAKT[42], and LTMTI[2].

####  3.3.3 Evaluation Metrics and Validation

We used the accuracy (ACC) and the area under the curve (AUC) as performance
metrics to compare the models’ performance in five datasets. We also used the
training speed, speed ratio, and memory usage as metrics to compare the
performance in the large-scale datasets containing long-sequence data (i.e.,
EdNet and Junyi Academy). Moreover, we used five-fold cross-validation for the
evaluation.

####  3.3.4 Hyperparameters for Experiments

To compare with each model, the same parameters were used for model training.
The batch size was set to 64, and the train/test split was 0.8/0.2. The model
used an embedding size of 128 and the Adam optimizer with a learning rate of
0.001. The loss function used was the Binary Cross Entropy with Logits Loss
(BCEWithLogitsLoss). The scheduler was set to OneCycleLR with a maximum
learning rate of 0.002. Dropout was also being used at a rate of 0.2. The
training ran for a total of 100 epochs, with early stopping set to 10 epochs.
If the validation loss does not decrease for the first three epochs, the
training stops, in order to prevent overfitting and save resources. The
maximum sequence length was 200, with an eight-attention head. Hidden sizes
were 128 for BERT, 512 for FFN, and 128 for LSTM. The Transformer
block/encoder layer was set to 12.

##  4 Results and Discussion

###  4.1 Overall Performance

LBKT outperforms four baseline models on most metrics in the experiments on
five benchmark datasets. Tabel 3 shows the overall performance of each model.
We used five-fold cross-validation to estimate their performances. LBKT
performed the best on EdNet and Junyi Academy datasets on both ACC and AUC
metrics. It also achieved the best performance on the ACC metric on assist12
and AUC on assist17. On algebra06, AKT achieved the best performance on the
ACC metric, BEKT achieved the best performance on the AUC metric, and LBKT
achieved the second-best performance on both metrics. This result indicates
that LBKT is an efficient KT model on most datasets, especially large-scale
datasets containing long-sequence interaction data. This was affected by our
LBKT model’s unique architecture. The LSTM block enables the model to learn
the sequential features of the long sequence and gives more importance to the
recent actions of the students, which prevents the model from giving too much
weight to the long-ago and low-relevance actions and thus improving the
training efficiency.

Table 3: Comparison of different KT models on five benchmark datasets. The best performance is denoted in bold. Dataset | Metrics | LBKT | BEKT | SSAKT | LTMTI | AKT | DKVMN
---|---|---|---|---|---|---|---
assist12 | ACC | 0.815 | 0.786 | 0.675 | 0.813 | 0.769 | 0.756
| AUC | 0.768 | 0.812 | 0.741 | 0.785 | 0.753 | 0.701
assist17 | ACC | 0.792 | 0.795 | 0.771 | 0.796 | 0.733 | 0.797
| AUC | 0.814 | 0.801 | 0.735 | 0.683 | 0.803 | 0.709
algebra06 | ACC | 0.801 | 0.797 | 0.795 | 0.811 | 0.831 | 0.800
| AUC | 0.799 | 0.815 | 0.774 | 0.791 | 0.814 | 0.793
EdNet | ACC | 0.803 | 0.781 | 0.761 | 0.799 | 0.756 | 0.800
| AUC | 0.815 | 0.795 | 0.798 | 0.802 | 0.798 | 0.796
Junyi | ACC | 0.832 | 0.807 | 0.777 | 0.797 | 0.791 | 0.790
Academy | AUC | 0.851 | 0.831 | 0.845 | 0.812 | 0.799 | 0.769

Tabel 6 shows the performance comparison on the two large-scale datasets. On
both datasets, LBKT achieved the best training efficiency. It was 4.29x faster
than BEKT on EdNet and 4.77x faster than BEKT on Junyi Academy. Compared with
the second-best model, AKT, LBKT was 1.32x faster on EdNet and 1.42x faster on
Junyi Academy. For the memory cost, LBKT was about one-third of BEKT and lower
than LTMTL on both datasets. Although the memory cost of LBKT was not the
smallest, LBKT has achieved the best results in both ACC and AUC metrics
running on the same GPU. This allows LBKT to run on middle-range GPUs. To
improve the training efficiency, we used a last input as the query method [5]
in the Transformer block instead of the whole sequence, which decreased the
complexity of the encoder to improve training speed and reduce memory cost.

Table 4: Performance comparison on the two large-scale datasets, EdNet and Junyi Academy. The best performance is denoted in bold. Model |  | EdNet |  |  | Junyi Academy |
---|---|---|---|---|---|---
| speed $\uparrow$ | speed ratio $\uparrow$ | memory $\downarrow$ | speed$\uparrow$ | speed ratio $\uparrow$ | memory $\downarrow$
BEKT | 4.93 | 1.00x | 16.7 GB | 4.85 | 1.00x | 16.6 GB
SSAKT | 7.13 | 1.44x | 3.4 GB | 6.22 | 1.28x | 3.2 GB
LTMTI | 13.8 | 1.32x | 7.69 GB | 12.1 | 1.19x | 8.82 GB
AKT | 17.1 | 3.25x | 4.32 GB | 16.4 | 3.35x | 4.37 GB
DKNMN | 5.97 | 2.34x | 7.68 GB | 4.67 | 3.75x | 8.53 GB
LBKT | 21.3 | 4.29x | 6.09 GB | 22.2 | 4.77x | 6.08 GB

To estimate the model performance on different lengths of data sequences, we
sorted the data in the two datasets by length and divided them into four sub-
datasets according to the average length. The average lengths of the four sub-
datasets were set as 100, 200, 300, and 400, respectively. Sequences shorter
than the mean were padded with 0s, and sequences longer than the mean were
pruned.

Fig. 2 shows the results of the speed performance comparison of each model
processing different lengths of data sequences. LBKT has a relatively high-
speed performance compared to other KT models when processing data sequences
with varying lengths. LBKT is the fastest model in all four groups of data
lengths. AKT and DKMN also have relatively high speeds, with AKT being the
second-fastest model in all groups and DKMN being the third-fastest model.
Overall, the results suggest that LBKT is the fastest model, and that it is
particularly efficient at dealing with long sequences of data. The fact that
LBKT maintains its high speed even when processing longer sequences of data
indicates that it is well-suited for tasks that require the analysis of large
amounts of data over extended periods of time.

Figure 2: Speed performance comparison of each model when processing data sequences with varying lengths. The vertical axis is the speed ($10^{4}$ samples per sec). Table 5: Results of the ablation study. LBKT-Rasch denotes LBKT without Rasch embedding; LBKT-LTSM denotes LBKT without LSTM block; and BERT denotes only the transformer structure-based blocks are included. The best performance is denoted in bold. Dataset | Metrics | LBKT | LBKT-Rasch | LBKT-LSTM | BERT
---|---|---|---|---|---
assist12 | ACC | 0.804 | 0.785 | 0.799 | 0.793
| AUC | 0.768 | 0.768 | 0.783 | 0.750
assist17 | ACC | 0.784 | 0.792 | 0.782 | 0.792
| AUC | 0.814 | 0.709 | 0.779 | 0.799
algebra06 | ACC | 0.801 | 0.796 | 0.792 | 0.798
| AUC | 0.799 | 0.756 | 0.809 | 0.765
EdNet | ACC | 0.803 | 0.729 | 0.722 | 0.801
| AUC | 0.815 | 0.758 | 0.794 | 0.809
Junyi Academy | ACC | 0.882 | 0.856 | 0.874 | 0.879
| AUC | 0.907 | 0.893 | 0.877 | 0.901

###  4.2 Ablation Study

In this section, we explore why LBKT performed better than other methods and
which components affected the overall performance. Table 5 shows the results
of the ablation study. We compared LBKT, LBKT without Rasch model-based
embeddings block (denoted as LBKT-Rasch), LBKT without LSTM block (denoted as
LBKT-LSTM), and LBKT without both Rasch model-based embeddings LSTM (denoted
as BERT). The results show that LBKT achieved the best performance on EdNet
and Junyi Academy on both ACC and AUC metrics. It also achieved the best
performance on one metric in every dataset. BERT-only achieved the best
performance on assist17 on ACC, which shows that the combination with Rasch
embeddings and LSTM could improve the performance of a single BERT model.

Table 4 shows the ablation study of speed performance comparison on the two
large-scale datasets. The results show that LBKT has the highest speed
performance among all models on both datasets, with a speed of 21.3 samples
per second on EdNet and 22.2 samples per second on Junyi Academy. This
suggests that LBKT is a highly efficient model for processing large amounts of
data in real-time. Interestingly, LBKT-LSTM, which removes the LSTM layer from
the proposed model, has a significantly lower speed performance compared to
LBKT, with a speed ratio of only 1.29x on EdNet and 1.19x on Junyi Academy.
This suggests that the LSTM layer is an important component in the proposed
model and contributes significantly to its speed performance. This is likely
due to the ability of LSTM to capture long-term dependencies and sequential
patterns in the data, which can be crucial in educational applications.

Table 6: Speed performance comparison ablation study on the two large-scale datasets, EdNet and Junyi Academy. The best performance is denoted in bold. Model |  | EdNet |  |  | Junyi Academy |
---|---|---|---|---|---|---
| speed | speed ratio | memory | speed | speed ratio | memory
LBKT | 21.3 | 4.29x | 6.09GB | 22.2 | 4.77x | 6.08GB
LBKT-Rash | 19.1 | 3.51x | 5.6GB | 20.31 | 3.67x | 3.1 GB
LBKT-LSTM | 13.8 | 1.29x | 7.87GB | 12.1 | 1.09x | 8.37GB
BERT | 12.1 | 3.27x | 4.32GB | 11.4 | 3.52x | 4.34GB

###  4.3 Analysis of Embedding Strategy

In this section, We used t-SNE as the visualisation tool to show the
interpretability of LBKT’s embedding strategy. Fig. 3-left shows the results
of No-Rasch-embedding, and Fig.3-right shows the Rasch embedding strategy. We
can see that, in the No-Rasch-embedding scenario, the difficult questions’
embeddings (dark blue vectors) mixed with the easy questions’ embeddings
(yellow to light blue vectors). In figure 3-right, the difficult level
embeddings were separated to avoid mixing with easy level embeddings.

Questions at a higher difficulty level are typically associated with longer
sequence data, as students spend more time and steps on difficult exercises,
which results in longer interaction sequences. Rasch model-based embeddings
could divide different difficulty-level parts before the start of the model
training and not mix them with other difficulty-level embeddings. As a result,
it might increase training efficiency to converge faster.

Figure 3: Visualisation of the embedding vector using t-SNE: without Rasch
embeddings (on the left) and with Rasch embeddings (on the right). The colour
bar is the predicted probability of the outputs.

##  5 Conclusion and Future work

In this study, we have developed LBKT, which employs a BERT-based architecture
with an LSTM block for processing long-sequence data, and Rasch model-based
embeddings for different difficulty levels of questions. Experiments show that
LBKT outperforms baseline models on most benchmark datasets. We also conducted
the speed performance experiment on the two large-scale datasets containing
long-sequence data. The results suggest that LBKT could process long-sequence
data faster and is more resource-efficient. Moreover, we conducted an ablation
study for different components of LBKT. The results indicate that the LSTM
component aided in improving the performance of dealing with long-sequence
data. Furthermore, we conducted an analysis of the embedding strategy using
t-SNE. The result shows that Rasch embedding could process the difficulty-
level features effectively. In future work, we plan to improve the
architecture to process more comprehensive data to satisfy the sustainable
development requirement of increasing large-scale datasets containing long-
sequence data emerging. Furthermore, we plan to build a more advanced
embedding strategy to allow the model to process multiple data types
simultaneously, such as exercises and concepts information from the dataset.

## Acknowledgments

This was was supported in part by……

## References

  * [1] Ghodai Abdelrahman, Qing Wang, and Bernardo Pereira Nunes.  Knowledge tracing: A survey.  ACM Computing Surveys, 2022.
  * [2] Youngduck Choi, Youngnam Lee, Dongmin Shin, Junghyun Cho, Seoyon Park, Seewoo Lee, Jineon Baek, Chan Bae, Byungsoo Kim, and Jaewe Heo.  Ednet: A large-scale hierarchical dataset in education.  In International Conference on Artificial Intelligence in Education, pages 69–73. Springer, 2020.
  * [3] Haw-Shiuan Chang, Hwai-Jung Hsu, and Kuan-Ta Chen.  Modeling exercise relationships in e-learning: A unified approach.  In EDM, pages 532–535, 2015.
  * [4] Yang Liu, Jing Zhou, and Weiguo Lin.  Efficient attentive knowledge tracing for long-tail distributed records.  In 2021 IEEE/ACIS 6th International Conference on Big Data, Cloud Computing, and Data Science (BCD), pages 104–109. IEEE, 2021.
  * [5] SeungKee Jeon.  Last query transformer rnn for knowledge tracing.  arXiv preprint arXiv:2102.05038, 2021.
  * [6] Qi Liu, Shuanghong Shen, Zhenya Huang, Enhong Chen, and Yonghe Zheng.  A survey of knowledge tracing.  arXiv preprint arXiv:2105.15106, 2021.
  * [7] Aritra Ghosh, Neil Heffernan, and Andrew S Lan.  Context-aware attentive knowledge tracing.  In Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pages 2330–2339, 2020.
  * [8] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  arXiv preprint arXiv:1810.04805, 2018.
  * [9] Georg Rasch.  Probabilistic models for some intelligence and attainment tests. ERIC, 1993.
  * [10] Ghodai Abdelrahman, Qing Wang, and Bernardo Nunes.  Knowledge tracing: A survey.  ACM Computing Surveys, 55(11):1–37, 2023.
  * [11] Albert T Corbett and John R Anderson.  Knowledge tracing: Modeling the acquisition of procedural knowledge.  User modeling and user-adapted interaction, 4(4):253–278, 1994.
  * [12] Albert Corbett.  Cognitive mastery learning in the act programming tutor.  In Adaptive User Interfaces. AAAI SS-00-01. Retrieved from https://aaai. org/Library/Symposia/Spring/ss00-01. php, 2000.
  * [13] Michael Villano.  Probabilistic student models: Bayesian belief networks and knowledge space theory.  In International Conference on Intelligent Tutoring Systems, pages 491–498. Springer, 1992.
  * [14] Michael V Yudelson, Kenneth R Koedinger, and Geoffrey J Gordon.  Individualized bayesian knowledge tracing models.  In International conference on artificial intelligence in education, pages 171–180. Springer, 2013.
  * [15] Tanja Käser, Severin Klingler, Alexander G Schwing, and Markus Gross.  Dynamic bayesian networks for student modeling.  IEEE Transactions on Learning Technologies, 10(4):450–462, 2017.
  * [16] Hao Cen, Kenneth Koedinger, and Brian Junker.  Learning factors analysis–a general method for cognitive model evaluation and improvement.  In Intelligent Tutoring Systems: 8th International Conference, ITS 2006, Jhongli, Taiwan, June 26-30, 2006. Proceedings 8, pages 164–175. Springer, 2006.
  * [17] Philip I Pavlik Jr, Hao Cen, and Kenneth R Koedinger.  Performance factors analysis–a new alternative to knowledge tracing.  Online Submission, 2009.
  * [18] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [19] Chris Piech, Jonathan Bassen, Jonathan Huang, Surya Ganguli, Mehran Sahami, Leonidas J Guibas, and Jascha Sohl-Dickstein.  Deep knowledge tracing.  Advances in neural information processing systems, 28, 2015.
  * [20] Shalini Pandey and George Karypis.  A self-attentive model for knowledge tracing.  arXiv preprint arXiv:1907.06837, 2019.
  * [21] Dongmin Shin, Yugeun Shim, Hangyeol Yu, Seewoo Lee, Byungsoo Kim, and Youngduck Choi.  Saint+: Integrating temporal features for ednet correctness prediction.  In LAK21: 11th International Learning Analytics and Knowledge Conference, pages 490–496, 2021.
  * [22] Hiromi Nakagawa, Yusuke Iwasawa, and Yutaka Matsuo.  Graph-based knowledge tracing: modeling student proficiency using graph neural network.  In IEEE/WIC/ACM International Conference on Web Intelligence, pages 156–163, 2019.
  * [23] Unggi Lee, Yonghyun Park, Yujin Kim, Seongyune Choi, and Hyeoncheol Kim.  Monacobert: Monotonic attention based convbert for knowledge tracing.  arXiv preprint arXiv:2208.12615, 2022.
  * [24] Zejie Tiana, Guangcong Zhengc, Brendan Flanaganb, Jiazhi Mic, and Hiroaki Ogatab.  Bekt: Deep knowledge tracing with bidirectional encoder representations from transformers.  In Proceedings of the 29th International Conference on Computers in Education, 2021.
  * [25] Katikapalli Subramanyam Kalyan, Ajit Rajasekharan, and Sivanesan Sangeetha.  Ammus: A survey of transformer-based pretrained models in natural language processing.  arXiv preprint arXiv:2108.05542, 2021.
  * [26] Niki Parmar, Ashish Vaswani, Jakob Uszkoreit, Lukasz Kaiser, Noam Shazeer, Alexander Ku, and Dustin Tran.  Image transformer.  In International Conference on Machine Learning, pages 4055–4064. PMLR, 2018.
  * [27] Luciano Floridi and Massimo Chiriatti.  Gpt-3: Its nature, scope, limits, and consequences.  Minds and Machines, 30(4):681–694, 2020.
  * [28] Zi-Hang Jiang, Weihao Yu, Daquan Zhou, Yunpeng Chen, Jiashi Feng, and Shuicheng Yan.  Convbert: Improving bert with span-based dynamic convolution.  Advances in Neural Information Processing Systems, 33:12837–12848, 2020.
  * [29] Fei Sun, Jun Liu, Jian Wu, Changhua Pei, Xiao Lin, Wenwu Ou, and Peng Jiang.  Bert4rec: Sequential recommendation with bidirectional encoder representations from transformer.  In Proceedings of the 28th ACM international conference on information and knowledge management, pages 1441–1450, 2019.
  * [30] Chris Donahue, Huanru Henry Mao, Yiting Ethan Li, Garrison W Cottrell, and Julian McAuley.  Lakhnes: Improving multi-instrumental music generation with cross-domain pre-training.  arXiv preprint arXiv:1907.04868, 2019.
  * [31] Zhaoxing Li, Mark Jacobsen, Lei Shi, Yunzhan Zhou, and Jindi Wang.  Broader and deeper: A multi-features with latent relations bert knowledge tracing model.  In European Conference on Technology Enhanced Learning, pages 183–197. Springer, 2023.
  * [32] Zhaoxing Li, Lei Shi, Alexandra Cristea, Yunzhan Zhou, Chenghao Xiao, and Ziqi Pan.  Simstu-transformer: A transformer-based approach to simulating student behaviour.  In International Conference on Artificial Intelligence in Education, pages 348–351. Springer, 2022.
  * [33] Zhaoxing Li, Lei Shi, Yunzhan Zhou, and Jindi Wang.  Towards student behaviour simulation: a decision transformer based approach.  In International Conference on Intelligent Tutoring Systems, pages 553–562. Springer, 2023.
  * [34] Zhaoxing Li.  Deep Reinforcement Learning Approaches for Technology Enhanced Learning.  PhD thesis, Durham University, 2023.
  * [35] Zhaoxing Li, Lei Shi, Jindi Wang, Alexandra I Cristea, and Yunzhan Zhou.  Sim-gail: A generative adversarial imitation learning approach of student modelling for intelligent tutoring systems.  Neural Computing and Applications, 35(34):24369–24388, 2023.
  * [36] Jindi Wang, Ioannis Ivrissimtzis, Zhaoxing Li, Yunzhan Zhou, and Lei Shi.  Exploring the potential of immersive virtual environments for learning american sign language.  In European Conference on Technology Enhanced Learning, pages 459–474. Springer, 2023.
  * [37] Jindi Wang, Ioannis Ivrissimtzis, Zhaoxing Li, and Lei Shi.  Comparative efficacy of 2d and 3d virtual reality games in american sign language learning.  In The 31st IEEE Conference on Virtual Reality and 3D User Interfaces. Newcastle University, 2024.
  * [38] Jindi Wang, Ioannis Ivrissimtzis, Zhaoxing Li, Yunzhan Zhou, and Lei Shi.  User-defined hand gesture interface to improve user experience of learning american sign language.  In International Conference on Intelligent Tutoring Systems, pages 479–490. Springer, 2023.
  * [39] Jindi Wang, Ioannis Ivrissimtzis, Zhaoxing Li, and Lei Shi.  Impact of personalised ai chat assistant on mediated human-human textual conversations: Exploring female-male differences.  In Companion Proceedings of the 29th International Conference on Intelligent User Interfaces, pages 78–83, 2024.
  * [40] Weicong Tan, Yuan Jin, Ming Liu, and He Zhang.  Bidkt: Deep knowledge tracing with bert.  In International Conference on Ad Hoc Networks, International Conference on Testbeds and Research Infrastructures, pages 260–278. Springer, 2022.
  * [41] Xia Sun, Xu Zhao, Yuan Ma, Xinrui Yuan, Feijuan He, and Jun Feng.  Muti-behavior features based knowledge tracking using decision tree improved dkvmn.  In Proceedings of the ACM Turing Celebration Conference-China, pages 1–6, 2019.
  * [42] Xuelong Zhang, Juntao Zhang, Nanzhou Lin, and Xiandi Yang.  Sequential self-attentive model for knowledge tracing.  In International Conference on Artificial Neural Networks, pages 318–330. Springer, 2021.