This draft was prepared using LaTeX template provided by Journal of Applied
Physics

# Droplet spreading: Role of characteristic length scales and modified
disjoining pressure

Debanik Bhattacharjee  Hadi Nazaripoor  Mohtada Sadrzadeh
<EMAIL_ADDRESS>Advanced Water Research Lab, Department of Mechanical
Engineering, University of Alberta, AB-T6G 2R3

###### Abstract

The numerical prediction of droplet base radius and contact angle depends on
the choice of characteristic radius and height of droplet. In this study, a
developed model on the basis of lubrication approximation is used to
investigate the effect of characteristic parameters on the spreading of
droplets over substrates. Based on two main stages of spreading (initial and
equilibrium), characteristic radius and height were first evaluated. The model
predictions were then compared with experimental results in literature for
both impermeable and permeable substrates. The study provides evidence that
the choice of characteristic length scales based on either stage enables
accurate prediction of the droplet base radius and contact angle. In addition,
a new scaling relation for the predictions of the numerical disjoining
pressure, with an error of $\pm$5%, is proposed.

##  I INTRODUCTION

Droplet spreading is ubiquitous in nature and has been exploited in many
applications such as microfluidic devices, Choi  _et al._ (2012) printing
processes, Kuang, Wang, and Song (2014) spray cooling, Mousa  _et al._ (2009)
and tissue assays. Jia and Qiu (2003) The spreading of a liquid on a substrate
is dictated by physical parameters,Benintendi and Smith (1999); Kovalchuk,
Trybala, and Starov (2014); Erbil and Meric (1997) topography, McHale, Newton,
and Shirtcliffe (2009); Courbin  _et al._ (2007); Savva and Kalliadasis (2012)
and chemical heterogeneities. Modaressi and Garnier (2002); Bliznyuk  _et al._
(2009); Herminghaus, Brinkmann, and Seemann (2008) The spreading behaviour of
a perfectly wetting liquid on a smooth, impermeable substrate is governed by
Tanner’s law that describes the time evolution of droplet radius ($\sim$
time1/10) and droplet height ($\sim$ time-1/5).Tanner (1979) Permeability of
the substrate is another pertinent feature which affects the spreading of the
droplet. Frank and Perre (2012); Acton, Huppert, and Worster (2001) Permeation
plays a significant role in the latter stages of spreading which is marked by
retraction of the contact line and consequently, complete imbibition of the
droplet volume. The characteristics of a perfectly wetting liquid spreading on
a dry, permeable substrate fall on two universal curves.Starov (2004) However,
no such universal behavior has been observed for partially wetting liquids on
permeable substrates.

Droplet spreading on a substrate is a classical example of a moving contact
line problem. Numerical techniques have been found as a versatile tool to
track the interface height in droplet spreading studies. Some of the most
commonly used techniques include Molecular Dynamics simulation,He and
Hadjiconstantinou (2003); De Ruijter, Blake, and De Coninck (1999); Song, Li,
and Liu (2013) Lattice Boltzmann method,Raiskinmäki  _et al._ (2000); Attar
and Körner (2009); Dupuis and Yeomans (2004) and lubrication based
approach.Ehrhard and Davis (1991); Espín and Kumar (2015); Benintendi and
Smith (1999); Zadražil, Stepanek, and Matar (2006) The scaling factors in the
lubrication based approach are radius and height of the droplet. The tri-phase
contact line is prone to shear-stress singularity due to the no-slip condition
which has been addressed through a number of approaches using lubrication
approximation. The most common approach is to use a slip model between the
liquid and the solid to relieve the stress singularity. Using this model
results in a nonlinear equation which can be solved either asymptotically
Greenspan (1978); Cox (1986); Haley and Miksis (1991); Bostwick (2013);
Ehrhard and Davis (1991) or numerically. Benintendi and Smith (1999); Smith
(1995) A constitutive relation describes the contact line velocity in terms of
the apparent contact angle. Although slip model provides accurate results, its
main shortcoming, however, lies in the accurate prediction of the
proportionality constant and the spreading exponent in the constitutive
relation since the spreading exponent varies for different liquid-substrate
combinations.McHale, Newton, and Shirtcliffe (2009)

Another approach to remove the contact line singularity is precursor film
assumption. The contact line region is described by disjoining pressure which
consists of attractive and repulsive components. Gomba and Homsy (2010);
Zadražil, Stepanek, and Matar (2006); Espín and Kumar (2015) The precursor
film or “primary film” was first observed for the water and acetic acid drops
spreading on glass and steel substrates.Hardy (1919) The thickness of the
precursor film was estimated to be few hundreds angstrom using ellipsiometry
and interferometry techniques, Bascom, Cottington, and Singleterry (1964)
scanning electron microscopy,Radigan  _et al._ (1974) and polarised reflection
microscopy.Ghiradella, Radigan, and Frisch (1975) The main challenge of this
approach is the sensitivity of the equilibrium contact angle to the disjoining
pressure parameter.

Identifying a proper characteristic radius and height is the first step
towards development of a predictive model for droplet spreading over
substrates. The characteristic values are typically selected at the onset of
equilibrium from experimental data.Espín and Kumar (2015) A possible
explanation for selecting the equilibrium characteristic values might be
related to the experimental measurement and easier evaluation of equilibrium
values. However, these studies would have been more interesting if the
sensitivity of the results to the selection of other characteristic values is
included. For a partially wetting liquid, the equilibrium contact angle has a
finite non-zero value and the equilibrium stage is interpreted to be
controlled by disjoining pressure. The disjoining parameter,
$\Pi=A\Big{[}(\frac{b}{h})^{n}-(\frac{b}{h})^{m}\Big{]}$ is a combination of
repulsive (first term) and attractive (second term) components acting when the
droplet height is small and comparable to the precursor film thickness, $b$,
$(n,m)=(3,2)$, and $A$ is the disjoining pressure parameter. It has been shown
that $A$ is inversely proportional to the precursor film thickness. It also
depends on the equilibrium contact angle and $(n,m)$.Schwartz (1998) However,
none of the proposed relations in the literature provide a viable $A$ which
leads to a better agreement between numerical prediction and experimentally
measured contact angle.

The present paper is aimed at investigating the effect of characteristic
length scales (radius and height) on the spreading of partially wetting
droplet over both impermeable and permeable substrates. A mathematical model
is used to describe the contact line region through precursor film
assumption.Espín and Kumar (2015) Using the no-slip condition at the
substrate, a nonlinear partial differential equation is derived which
describes the spatio-temporal evolution of droplet height based on lubrication
approximation. A disjoining pressure model is employed to remove contact line
singularity.Gomba and Homsy (2010); Zadražil, Stepanek, and Matar (2006);
Alleborn and Raszillier (2004) For imbibition, Poiseuille flow driven by
excess liquid pressure is assumed in each unconnected vertical pore which
essentially removes the need for thickness-dependent permeability.Davis and
Hocking (1999) Normalising factors are defined based on the initial and
equilibrium stages of droplet spreading. Using reported experimental data,De
Ruijter  _et al._ (1997); Kumar and Deshpande (2006); Wang  _et al._ (2016)
the aim is to investigate whether the prediction of droplet base radius and
contact angle is influenced by the stage-dependent characteristic radius and
height. In addition, the numerical disjoining pressure parameter is predicted
based on its theoretical value through different stages of spreading.

Section IIA includes details about the governing equations used. Section IIB
presents boundary conditions needed for solving the resulting partial
differential equation. In section IIC, the variables which are used to non-
dimensionalise the governing equations and boundary conditions, are presented.
The derivation of the time evolution equation is discussed in section IID
along with the numerical scheme utilised to solve the equation. The
normalising factors are defined in section IIE. In section III, a comparison
of the numerical predictions with the experimental data is presented.De
Ruijter  _et al._ (1997); Kumar and Deshpande (2006); Wang  _et al._ (2016) In
section IV, a new scaling to bridge the gap between the theoretical and model
value of the disjoining pressure parameter is proposed for both impermeable
and permeable substrates. All findings are summarised in section V.

##  II MATHEMATICAL FORMULATION

Figure 1: (a) Geometrical configuration of an axi-symmetric droplet and
different stages of droplet spreading on (b) an impermeable substrate and (c)
a permeable substrate. The subscript ’in’ and ’eq’ denote the initial and
equilibrium value of a quantity, respectively.

An axi-symmetric liquid droplet is considered to spread on a smooth,
horizontal substrate (Fig. 1(a)). The dynamics of droplet spreading depends on
the nature of the substrate and is different for impermeable and permeable
substrate, as illustrated in Fig. 1(b) and Fig. 1(c). The droplet is assumed
to be an incompressible, Newtonian liquid with constant density $\rho$,
dynamic viscosity, $\mu$ and surface tension, $\sigma$. Since the droplet is
assumed to be axi-symmetric in nature, the governing equations are considered
in cylindrical coordinates.

###  II.1 GOVERNING EQUATIONS

To describe the droplet spreading, pressure, $p$ and velocity field,
$\textbf{v}=(u,w)$ is introduced where $u$ and $w$ denote the velocity in
radial and vertical directions respectively. The continuity equation and
Navier-Stokes equations are as follows:

$\nabla.\textbf{v}=0\\\ $ (1) $\rho(\frac{\partial\textbf{v}}{\partial
t}+\textbf{v}.\nabla\textbf{v})=\mu\Delta\textbf{v}-\nabla p-\rho\textbf{g}\\\
$ (2)

where g is the acceleration of gravity.

###  II.2 BOUNDARY CONDITIONS

At the liquid-gas interface, $z=h(r,t)$, the kinematic condition is used to
relate the velocity field to interface height.

$w=\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial r}$ (3)

The normal and tangential component of interfacial stress balance is as
follows:

$\hat{\textbf{n}}.\textbf{T}.\hat{\textbf{n}}=2H\sigma-\Pi,\;\hat{\textbf{t}}.\textbf{T}.\hat{\textbf{n}}=t.\nabla\sigma$
(4)

where $\textbf{T}=\mu(\nabla u+\nabla u^{T})$ is the viscous stress of the
liquid, $\hat{\textbf{n}}$ and $\hat{\textbf{t}}$ are unit vector normal and
tangential to the liquid-gas interface respectively which are defined as
follows:

$\hat{\textbf{n}}=\frac{(-\frac{\partial h}{\partial
r},1)}{\sqrt{1+(\frac{\partial h}{\partial
r})^{2}}},\hat{\textbf{t}}=\frac{(1,\frac{\partial h}{\partial
r})}{\sqrt{1+(\frac{\partial h}{\partial r})^{2}}}\\\ $ (5)

The curvature of the substrate is given by

$2H=-({\frac{\partial^{2}h}{\partial r^{2}}}+\frac{1}{r}\frac{\partial
h}{\partial r})\\\ $ (6)

A two term disjoining pressure, $\Pi$ Alleborn and Raszillier (2004);
Schwartz, Roux, and Cooper-White (2005) which is used to relieve the stress
singularity at the contact line is given by

$\Pi=A\Big{[}(\frac{b}{h})^{n}-(\frac{b}{h})^{m}\Big{]}\\\ $ (7)

where $A\geq 0$ is the disjoining pressure parameter and exponents,
$(n,m)=(3,2)$ are constant and $b$ is the precursor film thickness. The value
of $A$ is related to the equilibrium contact angle, $\theta_{eq}$ as
follows:Schwartz, Roux, and Cooper-White (2005)

$A=\frac{(n-1)(m-1)}{b(n-m)}\sigma(1-cos\theta_{eq})\\\ $ (8)

At the liquid-substrate interface, $z=0$, the imbibition condition is used:

$w=\frac{\kappa}{\mu}p\\\ $ (9)

where $\kappa$ is the proportionality constant with units of length and
depends on pore diameter, number of pores per unit area, and substrate
thickness. Davis and Hocking (1999) No-slip condition is as follows:

$u=0\\\ $ (10)

The following conditions are imposed for axi-symmetric spreading:

$\frac{\partial h}{\partial r}\bigg{|}_{r=0}=0,\frac{\partial^{3}h}{\partial
r^{3}}\bigg{|}_{r=0}=0$ (11)

For droplet volume $V_{0}$ to be conserved, following equation is used:

$\int\limits_{0}^{r_{0}}2\pi rh(r)dr=V_{0}+\pi b$ (12)

where subscript 0 denotes the initial value of a quantity and $\pi b$ is the
correction factor for precursor film on a circular region with radius equal to
the characteristic value.

###  II.3 SCALINGS

The governing equations and the corresponding boundary conditions are
normalized using following factors:

$r=l\tilde{r},z=h\tilde{z},u=u_{sc}\tilde{u},w=\epsilon u_{sc}\tilde{w}$
$t=\frac{l}{u_{sc}}\tilde{t},p=p_{sc}\tilde{p},\Pi=p_{sc}\tilde{\Pi}$ (13)

where $u_{sc}=\frac{\epsilon^{3}\sigma}{3\mu}$ is the characteristic spreading
speed, $p_{sc}=\frac{h\sigma}{l^{2}}$ is the pressure scale, $b=h\tilde{b}$ is
the precursor film thickness, and $\epsilon=\frac{h}{l}$ is the lubrication
ratio. The effective permeability of the substrate is given by
$\kappa=\frac{dh^{3}}{3l^{2}}\tilde{\kappa}$ where $d$ is the thickness of the
substrate. The non-dimensional quantities are denoted with tilde. The droplet
is assumed to be thin, which means the ratio of characteristic height, $h$, to
characteristic radius, $l$, is much less than unity ($\epsilon\ll 1$). Hence,
the lubrication approximation is applied which reduces the governing equations
and boundary conditions to:
Continuity equation:

$\frac{\partial{\tilde{w}}}{\partial{\tilde{z}}}+\frac{1}{\tilde{r}}\frac{\partial({\tilde{r}\tilde{u}})}{\partial\tilde{r}}=0,\\\
$ (14)

r-momentum equation:

$\frac{\partial^{2}{\tilde{u}}}{\partial{\tilde{z}}^{2}}-3\frac{\partial{\tilde{p}}}{\partial{\tilde{r}}}=0,$
(15)

z-momentum equation:

$-\frac{\partial{\tilde{p}}}{\partial{\tilde{z}}}-B=0$ (16)

where $B=\frac{l^{2}\rho g}{\sigma}$ is the Bond number.
The kinematic condition at $\tilde{z}=\tilde{h}(\tilde{r},\tilde{t})$, is as
follows:

$\tilde{w}=\frac{\partial\tilde{h}}{\partial\tilde{t}}+u\frac{\partial\tilde{h}}{\partial\tilde{r}},$
(17)

The normal and tangential component of interfacial stress balance is written
as follows:

$\tilde{p}=-(\frac{\partial^{2}{\tilde{h}}}{\partial{\tilde{r}^{2}}}+\frac{1}{\tilde{r}}\frac{\partial{\tilde{h}}}{\partial{\tilde{r}}})-\tilde{A}\Big{[}(\frac{\tilde{b}}{\tilde{h}})^{3}-(\frac{\tilde{b}}{\tilde{h}})^{2}\Big{]},$
(18) $\frac{\partial{\tilde{u}}}{\partial{\tilde{z}}}=0\\\ $ (19)

where $\tilde{A}=\frac{A}{p_{sc}}$.
Imbibition condition at $\tilde{z}=\tilde{b}$ is as follows:

$\tilde{w}=-\frac{\tilde{\kappa}}{\mu d}\tilde{p}\\\ $ (20)

No-slip condition at $\tilde{z}=0$ is as follows:

$\tilde{u}=0$ (21)

###  II.4 NUMERICAL SCHEME

The non-dimensional time-evolution equation for the droplet height is obtained
by integrating eqns.(14-16) subject to the boundary conditions (17-21):

$\frac{\partial{\tilde{h}}}{\partial{\tilde{t}}}=\frac{1}{\tilde{r}}\frac{\partial}{\partial{\tilde{r}}}\Bigg{\\{}\tilde{r}\tilde{h}^{3}\frac{\partial}{\partial{\tilde{r}}}\Bigg{(}B\tilde{h}-\tilde{A}{\tilde{\Pi}}-\frac{\partial^{2}\tilde{h}}{\partial\tilde{r}^{2}}-\frac{1}{\tilde{r}}\frac{\partial{\tilde{h}}}{\partial{\tilde{r}}}\Bigg{)}\Bigg{\\}}-\tilde{\kappa}\Bigg{\\{}B(\tilde{h}-\tilde{b})-\tilde{A}{\tilde{\Pi}}-\frac{\partial^{2}\tilde{h}}{\partial\tilde{r}^{2}}-\frac{1}{\tilde{r}}\frac{\partial{\tilde{h}}}{\partial{\tilde{r}}}\Bigg{\\}}$
(22)

The solution of eqn.(22) is subject to the following boundary conditions on
the domain $0\leq\tilde{r}\leq\textit{L}$:

$\frac{\partial{\tilde{h}}}{\partial{\tilde{r}}}(0,\tilde{t})=0,\tilde{h}(L,\tilde{t})=\tilde{b},\frac{\partial^{3}h}{\partial{\tilde{r}^{3}}}(0,\tilde{t})=0,\frac{\partial{\tilde{h}}}{\partial{\tilde{r}}}(L,\tilde{t})=0.$
(23)

The initial droplet profile is obtained in the form of a fourth-degree
polynomial satisfying the above boundary conditions (23) and the volume
constraint (12). Computational domain length, L is considered to be 4 with
600-700 points per unit length. The spatial derivatives are discretised using
second-order central finite difference scheme which converts the PDE into sets
of differential algebraic equations (DAE). For the resulting system of DAEs,
an adaptive time-step solver, DASSL, is used which is available in SLATEC
library.

The apparent contact angle is calculated by finding maximum slope of the
interface and multiplying by $\epsilon$ to convert into dimensional form. The
radius is obtained by finding the distance between the z-axis and the
intersection point of the apparent contact line with the r-axis.

###  II.5 NORMALISING FACTORS

In order to obtain eqn.(22), two normalising factors are used. The first one
is based on the radius and height of the droplet at the onset of equilibrium
and the second one is described by those values at the initial stage. Based on
this definition, the required parameters are calculated using both normalising
factors as described below:
a. First normalising factor: Using equilibrium values ($l_{eq}$,$h_{eq}$) as
normalising factors, the following input parameters are obtained: Normalised
initial radius, $r_{0}=\frac{l_{in}}{l_{eq}}$; normalised initial height,
$h_{0}=\frac{h_{in}}{h_{eq}}$; initial volume based on spherical-cap
approximation, $V_{in}=\frac{\pi h_{in}(3l_{in}^{2}+h_{in}^{2})}{6}$; final
volume when the droplet ceases to exist on the substrate (assuming spherical-
cap), $V_{eq}=\frac{\pi h_{eq}(3l_{eq}^{2}+h_{eq}^{2})}{6}$; normalised
initial volume, $V_{0}=\frac{V_{in}}{V_{eq}}$; Bond number,
$B_{eq}=\frac{l_{eq}^{2}\rho g}{\sigma}$; time
$t_{eq}=\frac{h_{eq}^{3}\sigma}{3l_{eq}^{4}\mu}\tilde{t}$; and lubrication
ratio, $\epsilon_{eq}=\frac{h_{eq}}{l_{eq}}$.
b. Second normalising factor: Using initial values ($l_{in}$,$h_{in}$) as
normalising factors, the following input parameters are found as:
$r_{0}=\frac{l_{in}}{l_{in}}=1$; $h_{0}=\frac{h_{in}}{h_{in}}=1$;
$V_{0}=\frac{V_{in}}{V_{in}}=1$; $B_{in}=\frac{l_{in}^{2}\rho g}{\sigma}$;
$t_{in}=\frac{h_{in}^{3}\sigma}{3l_{in}^{4}\mu}\tilde{t}$; and
$\epsilon_{in}=\frac{h_{in}}{l_{in}}$

##  III RESULTS AND DISCUSSION

For a small droplet, where the effect of gravity is negligible, the driving
force for spreading is capillary pressure. The maximum extent of spreading is
dictated by the dominant effect of disjoining pressure over capillary pressure
in the later stages. For a smooth, impermeable substrate, a perfectly wetting
Newtonian liquid obeys Tanner’s law.Tanner (1979) However, the maximum
spreading radius will vary if the liquid partially wets the substrate. Under
the combined effect of attractive and repulsive components of intermolecular
forces in disjoining pressure, a liquid partially wets the substrate and thus,
has a finite non-zero equilibrium contact angle.Glasner (2003) For smooth
permeable substrates, the spreading occurs at two stages: First, the inertial
stage where lubrication ratio is much larger than unity. At this stage, the
droplet spreads to a maximum radius and stops momentarily. Second, the
equilibrium stage which is signified by a decrease in droplet volume due to
the imbibition. As a result, the droplet radius decreases and the contact
angle remains constant over a particular time before becoming zero. At the
equilibrium stage, the lubrication ratio is much less than unity.
A 3-step procedure is followed in order to make comparison of model
predictions with experimental data:
(i) A specific value of $b$ which yields $\tilde{b}$ through the scaling is
used as a starting point.
(ii) Normalised initial radius ($r_{0}$) and height ($h_{0}$) are obtained
using both normalising factors. These parameters along with $\tilde{b}$ are
utilised to obtain the initial droplet profile for each case. The model
parameters used for matching with the experimental values are disjoining
pressure parameter, $\tilde{A}$ and non-dimensional imbibition coefficient,
$\tilde{k}$. $\tilde{A}$ governs the equilibrium stage of spreading and its
higher values signifies less hydrophilicity or higher equilibrium contact
angle. The equilibrium stage is first matched by varying $\tilde{A}$, then
$\tilde{k}$ is adjusted so as to match the imbibition time. Many combinations
of ($\tilde{A}$, $\tilde{k}$) are tried for a specific $\tilde{b}$ until good
agreement with the experimental results is found.
(iii) For both impermeable and permeable substrates, ($A,k$) corresponding to
each value of $b$ are obtained using the first normalising factor. The
simulation is then conducted using these values and the second normalising
factor. This approach provides a clear picture on the role of each normalising
factor on the numerical predictions.

Figure 2: Comparison between experimental data and simulation results of
contact angle of squalene spreading on glass.De Ruijter  _et al._ (1997) The
simulation parameters used are: $V_{in}=V_{eq}=6.046\;\mu$L, $l_{in}=1.6$ mm,
$l_{eq}=2.548$ mm, $\epsilon_{in}=0.781$, $\epsilon_{eq}=0.229$,
$B_{in}=0.702$, $B_{eq}=1.78$, $b=0.98\;\mu$m, $A=8271\;$N$m^{-2}$

In Fig. 2, the contact angle of squalene drop spreading on glass De Ruijter
_et al._ (1997) is compared with the model predictions. There is no change in
volume over time since glass is impermeable. The spherical-cap approximation
is used for calculating the volume of the droplet. First the initial values
are normalised with the equilibrium characteristic values. For a specific $b$
($<$1 $\mu$m Popescu  _et al._ (2012)), the initial shape of the droplet is
obtained. As can be seen in Fig. 2, there is a good agreement between model
predictions and experimental data, only at the equilibrium stage of spreading
and the model failed in predicting the earlier stages of spreading. Then,
using the corresponding ($b,A$) obtained by the first normalising factor, the
simulation is performed based on the second normalising factor. A clear
overlap in the numerical prediction highlights the fact that either
normalising factor can predict the equilibrium contact angle. The model was
further tested for $b=0.75\;\mu$m and $0.5\;\mu$m and a good match was
observed (The results are not provided since the behaviour was exactly the
same). For each $b$, the numerical predictions coincided using both
normalising factors.

Figure 3: (a) The nature of defects considered. Comparison between
experimental data and simulation results of (b) droplet base radius, and (c)
contact angle of squalene spreading on glass in the presence of substrate
roughness. De Ruijter  _et al._ (1997) The simulation parameters remain
unchanged.

The decrease in contact angle is sharper in numerical simulation compared to
the experiments and equilibrium is reached earlier. The difference between
numerically predicted results and experiments might be attributed to the
substrate roughness. Since any substrate possesses inherent roughness (unless
otherwise specified), one way to delay the onset of equilibrium is to take
into account the effect of substrate roughness in the mathematical model. In
the presence of roughness, the contact line is pinned and the apparent contact
angle increases.Lafuma and Quéré (2003) To incorporate roughness in the model,
a shape function, $\tilde{z}=\gamma(\tilde{r})$ is introduced in this study.
The governing equations remain unchanged in the presence of substrate
roughness. However, the kinematic condition and the interfacial stress balance
components are evaluated at
$\tilde{z}=\gamma(\tilde{r})+\tilde{h}(\tilde{r},\tilde{t})$. The imbibition
and no-slip conditions also occur at $\tilde{z}=\gamma(\tilde{r})+\tilde{b}$
and at $\tilde{z}=\gamma(\tilde{r})$ respectively. Based on these
modifications, the droplet evolution equation takes the following form:

$\displaystyle\frac{\partial{\tilde{h}}}{\partial{\tilde{t}}}=\frac{1}{\tilde{r}}\frac{\partial}{\partial{\tilde{r}}}\Bigg{\\{}\tilde{r}\tilde{h}^{3}\frac{\partial}{\partial{\tilde{r}}}\Bigg{(}B(\tilde{h}+\gamma)-\tilde{A}{\tilde{\Pi}}-\frac{\partial^{2}(\tilde{h}+\gamma)}{\partial\tilde{r}^{2}}-\frac{1}{\tilde{r}}\frac{\partial{(\tilde{h}+\gamma)}}{\partial{\tilde{r}}}\Bigg{)}\Bigg{\\}}$
$\displaystyle-\tilde{\kappa}\Bigg{\\{}B(\tilde{h}-\tilde{b})-\tilde{A}{\tilde{\Pi}}-\frac{\partial^{2}(\tilde{h}+\gamma)}{\partial\tilde{r}^{2}}-\frac{1}{\tilde{r}}\frac{\partial{(\tilde{h}+\gamma)}}{\partial{\tilde{r}}}\Bigg{\\}}$
(24)

There are complicated forms of shape function, $\gamma(r)$, in the literature
that are able to imitate a real rough substrate very accurately.Chow (1998);
Savva, Kalliadasis, and Pavliotis (2010) However, in this study, the most
simplified form of this function, known as Gaussian model,Herminghaus (2012)
is utilised. The Gaussian model of roughness is represented by the following
equation:

$\gamma(r)=-r_{1}exp\Bigg{(}-\frac{(r-r_{loc})^{2}}{2r_{2}^{2}}\Bigg{)}$ (25)

where $r_{1},r_{2}=O(\tilde{b})$ and $r_{loc}$ represents the location of the
defect.

For the present simulation, $r_{1}$ and $r_{2}$ are selected as
$50~{}\tilde{b}$ and $5~{}\tilde{b}$, respectively. To investigate the role of
defects on contact line motion, three different configurations are selected,
as shown in Fig. 3(a). When the defect is placed very close to the initial
shape (Defect (i) in Fig. 3), the pinning phenomenon is observed which is
represented by no change in radius over a certain time period. If the nature
of the defect remains unchanged but is placed much farther from the initial
shape (Defect (ii)), pinning persists for comparatively longer period. This is
attributed to non-uniform speed of droplet spreading. The characteristic speed
of a spreading droplet is much higher at the initial stages and as time passes
the spreading front slows down owing to substrate irregularities. Since the
height of the droplet also decreases over time, the rate of conversion of the
potential energy (owing to greater height in the initial stages) to kinetic
energy diminishes which is manifested in the speed of the spreading front. As
a result, for a defect placed much farther on the path, the contact line has a
tendency to stay pinned for a longer time and thus, the equilibrium is
delayed. If there is another identical defect placed at a certain distance
from the first defect (Defect (iii)), a staircase-like behaviour is observed
marked by a delayed equilibrium. Our preliminary investigations showed that
this staircase behaviour depends on the nature, density and spacing of
defects. The detailed study of staircase bahviour with the aim of acquiring
more realistic replication of a rough substrate is the topic of our ongoing
research.

Figure 4: Comparison between experimental data and simulation results of
contact angle of DBP spreading on PET.De Ruijter  _et al._ (1997) The
simulation parameters used are: $V_{in}=V_{eq}=5\;\mu$L, $l_{in}=1.454$ mm,
$l_{eq}=2.238$ mm, $\epsilon_{in}=0.839$, $\epsilon_{eq}=0.277$,
$B_{in}=0.635$, $B_{eq}=1.504$, $b=0.98\;\mu$m, $A=13587.2\;$N$m^{-2}$

In another study of impermeable substrate, the contact angle of di-n-butyl
phthalate (DBP) spreading on poly(ethyleneterephthalate) (PET) is
considered.De Ruijter  _et al._ (1997) A similar procedure as the previous
case is followed and a good agreement with the experimental data is achieved
for $b=0.98\;\mu$m (Fig. 4), $0.75\;\mu$m, and $0.5\;\mu$m. Comparing the
model predicted results and experiments show that irrespective to the applied
characteristic lengths (for radius and height) and precursor film, $b$, the
equilibrium contact angle is mainly governed by the proper choice of $A$.

Figure 5: Comparison between experimental data and simulation results of (a)
droplet base radius, and (b) contact angle of water droplet spreading on a
h-PBT-4 mat.Wang  _et al._ (2016) The simulation parameters used are:
$V_{in}=10.04\;\mu$L, $V_{eq}=1.089\;\mu$L, $l_{eq}=2.81$ mm, $l_{in}=2.087$
mm, $\epsilon_{in}=0.281$, $\epsilon_{eq}=0.076$, $B_{in}=1.076$,
$B_{eq}=0.594$, $b=0.98\;\mu$m, $A=6379.2\;$N$m^{-2}$, $k/d=5.076\times
10^{-8}$

For a permeable substrate, model predictions and experimental measurements of
water droplet spreading on a poly(butylene ter-ephthalate) (h-PBT-4) mat are
compared in the absence of contact line pinning for $b=0.98\;\mu$m (Fig. 5).
In the experimental work, the onset of lubrication phase was reported as
$t=0.117\;$s. The corresponding initial height ($h_{in}$) and radius
($l_{in}$) of the droplet were $0.79\;$mm and $2.81\;$mm, respectively. The
equilibrium height ($h_{eq}$) and radius ($l_{eq}$) were reported to be
$0.159\;$mm and $2.087\;$mm, respectively. The initial droplet profile is
obtained based on these values. As can be observed in Fig. 5, there is an
overlap between numerical predictions using either normalising factor for
constant ($b,A,k$). In both cases, droplet radius increases at the initial
stages followed by a retraction due to permeation which is similar to what is
observed in experiments. Using other $b$ values ($b=0.75\;\mu$m and
$0.5\;\mu$m), the same results are obtained. This reinforces our hypothesis
that using either normalising factors enables accurate prediction of the
experimental values, for both permeable and non-permeable substrates.
The time lag between the numerical prediction and experimental values is
primarily related to the starting point of simulation at $t=0.117\;$s as this
marks the onset of lubrication. A possible explanation might be the difference
between experimental and numerical initial droplet profile to run the
simulation.

To showcase another example of the spreading phenomena on a permeable
substrate, the droplet radius and contact angle of silicone oil grade SO300
spreading on unidirectional mat (UDMAT) Kumar and Deshpande (2006) are
compared with the model predictions. The initial and equilibrium height and
radius of the droplet were obtained at $t=0$ and at the onset of equilibrium
stage respectively. Using first normalising factor resulted in better
agreement with experimental data for $b=0.98\;\mu$m (Fig. 6) and $0.75\;\mu$m.
When compared to the second normalising factor, there is an offset in
prediction which might be due to the droplet profile employed in this case for
running simulation.

Figure 6: Comparison between experimental data and simulation results of (a)
droplet base radius, and (b) contact angle of SO300 droplet spreading on
UDMAT.Kumar and Deshpande (2006) The simulation parameters used are:
$V_{in}=7.788\;\mu$L, $V_{eq}=1.609\;\mu$L, $l_{in}=1.7$ mm, $l_{eq}=2.32$ mm,
$\epsilon_{in}=0.82$, $\epsilon_{eq}=0.082$, $B_{in}=1.404$, $B_{eq}=2.615$,
$b=0.98\;\mu$m, $A=9633$ N$m^{-2}$, $k/d=2.039\times 10^{-6}$

The significant difference between model prediction of the contact angle (Fig.
6(b)) for $t<2\;$s is attributed to large lubrication ratio ($\epsilon\gg 1$)
at the inertial stage. The second normalising factor hinges on inertial values
which was not measured accurately. Even though dynamic radius and height
measurements are exclusively reported, the literature did not provide any
proper explanation for the calculation of the contact angle.Kumar and
Deshpande (2006)

##  IV MODIFIED SCALING FOR BETTER PREDICTION OF EQUILIBRIUM CONTACT ANGLE

Theoretically, the disjoining pressure parameter, $A_{th}$ is related to the
equilibrium contact angle, $\theta_{eq}$ as follows:Schwartz, Roux, and
Cooper-White (2005)

$A_{th}=\frac{(n-1)(m-1)}{b(n-m)}\sigma(1-cos\theta_{eq})\;or,\frac{2}{b}\sigma(1-cos\theta_{eq})\\\
$ (26)

To numerically model the spreading process, it is essential to start with an
initial guess of $\tilde{A}$. Since the equilibrium contact angle,
$\theta_{eq}$, is known from the dynamic contact angle measurement, it is
possible to have an estimate of $\tilde{A}$ ($A=p_{sc}\tilde{A}$).

It is found that the suggested $A$ (from eqn.26) cannot predict the
equilibrium contact angle accurately. Hence, the following modification is
proposed:

$A_{m}=\frac{A_{th}}{\epsilon^{n}}$ (27)

where $n$ is the exponent to be ascertained for impermeable and permeable
substrates.

The role of both normalising factors on the prediction of $n$ for droplet
spreading on impermeable and permeable substrate is presented next.

a. First normalising factor:

(i) Impermeable cases: The squalene droplet spreading on glass is considered
as an impermeable case. Starting with the first normalising factor, $n$ is
varied from 0 to 1.15 until the ratio $A_{m}/A$ was close to 1, with an error
of $\pm$5$\%$. According to Fig. 7(a), choosing $n=1.12$ resulted in a better
starting guess of $\tilde{A}$ for a particular precursor film thickness.
Similarly, for DBP spreading on PET, $n$ is changed from 0 to 1.3 (Fig. 8(a)).
A choice of $n=1.24$ has led to a better starting guess of $\tilde{A}$ for
each value of precursor film thickness.

(ii) Permeable cases: For water droplet spreading on h-PBT-4 mat, $n$ is
obtained to be $0.95$, $0.8$, $0.7$ for $b=0.98\;\mu$m, $0.75\;\mu$m, and
$0.5\;\mu$m, respectively (Fig. 9(a)). Following the similar procedure for
SO300 droplet spreading on UDMAT, the value of $n$ is evaluated to be $0.65$
with an error within $\pm$5$\%$ for $b=0.98\;\mu$m and $0.75\;\mu$m (Fig.
10(a)).

Figure 7: Ratio $A_{m}/A$ for (a) First normalising factor and (b) Second
normalising factor of squalene spreading on glass with
$\theta_{eq}=12.88^{\circ}$.De Ruijter  _et al._ (1997) Figure 8: Ratio
$A_{m}/A$ for (a) First normalising factor and (b) Second normalising factor
of DBP spreading on PET with $\theta_{eq}=15.569^{\circ}$.De Ruijter  _et al._
(1997)

b. Second normalising factor:

(i) Impermeable cases: Figs. 7(b) and 8(b) show that, using second normalising
factor the value of $n$ are obtained to be $6.6$ and $9.2$ for squalene
droplet spreading on glass and DBP spreading on PET, respectively. It is
noteworthy that in both experiments, no effort was made to precisely track the
onset of lubrication. This could be one plausible explanation for such a big
variation in the value of $n$ for these impermeable cases.

Figure 9: Ratio $A_{m}/A$ for (a) First normalising factor and (b) Second
normalising factor of water droplet spreading on a h-PBT4 mat with
$\theta_{eq}=5^{\circ}$. Wang  _et al._ (2016) Figure 10: Ratio $A_{m}/A$ for
(a) First normalising factor and (b) Second normalising factor of SO300
spreading on UDMAT with $\theta_{eq}=17^{\circ}$.Kumar and Deshpande (2006)

(ii) Permeable cases: For water droplet spreading on h-PBT-4 mat, the values
of $n=1.9$, $1.6$, and $1.4$ are obtained for $b=0.98\;\mu$m, $0.75\;\mu$m,
and $0.5\;\mu$m, respectively (Fig. 9(b)). The onset of lubrication was
precisely tracked. The second normalising factor yields value of $n=8.5$ for
$b=0.98$ and $0.75\;\mu$m for SO300 spreading on UDMAT.

In summary, first normalising factor yields $n$ between 1.0-1.3 for
impermeable and 0.5-1.0 for permeable substrates, thereby leading to a better
starting guess for $\tilde{A}$. For second normalising factor, it is essential
to precisely track the initial conditions. Only one experimental study on a
permeable substrate meets this criterion.Wang  _et al._ (2016) Hence, in the
absence of profound evidence, it is difficult to specify a range of $n$ to be
used for the second normalising factor.

##  V CONCLUSION

Ability of a mathematical model to accurately predict the droplet radius and
apparent contact angle is governed by the choice of radius and height at the
initial and equilibrium stages of spreading. Using four examples of droplet
spreading over impermeable and permeable substrates, it is observed that
choosing radius and height in either stage of the spreading process provides
similar prediction of the experimentally obtained base radius and equilibrium
contact angles. In addition, incorporating the substrate roughness into the
model is found to bridge the gap between model predictions and experiments for
initial stages of spreading. A modified scaling which relates the theoretical
value of the disjoining pressure parameter to its numerical counterpart
through the lubrication ratio is also proposed. The lubrication ratio takes
different exponents based on the permeability of the substrate. For
impermeable and permeable substrates, $n$ was found to be in the range of
1.0-1.3 and 0.5-1.0 respectively. The proposed method greatly simplified the
initial guess for the disjoining pressure parameter in the numerical
simulation as previously there was no possibility of ascertaining the value.
This study has shown that it is always possible to predict the value of the
lubrication exponent using the first normalising factor which depends on the
equilibrium condition. However, the second normalising factor depends on the
initial conditions entirely and thus, without precise tracking of the
lubrication ratio it is not feasible to suggest a range of $n$ for
impermeable/permeable substrates. Overall, the equilibrium stage
characteristic parameters is more appropriate since the modified scaling
allows prediction of the numerical disjoining parameter with $\pm 5\%$ error.
As a future study, more effort is needed for the accurate calculation of the
initial droplet profile and shape approximation for volume at each stage of
spreading. The hypotheses provided in this work is expected to hold true when
the model is extended to include the whole spectrum of spreading, especially
the inertial stages when the lubrication ratio is much greater than unity.

## References

  * Choi  _et al._ (2012) K. Choi, A. H. Ng, R. Fobel, and A. R. Wheeler, Annu. Rev. Anal. Chem. 5, 413 (2012).
  * Kuang, Wang, and Song (2014) M. Kuang, L. Wang, and Y. Song, Adv. Mater. 26, 6950 (2014).
  * Mousa  _et al._ (2009) N. A. Mousa, M. J. Jebrail, H. Yang, M. Abdelgawad, P. Metalnikov, J. Chen, A. R. Wheeler, and R. F. Casper, Sci. Transl. Med. 1, 1ra2 (2009).
  * Jia and Qiu (2003) W. Jia and H.-H. Qiu, Exp. Therm Fluid Sci. 27, 829 (2003).
  * Benintendi and Smith (1999) S. W. Benintendi and M. K. Smith, Phys. Fluids 11, 982 (1999).
  * Kovalchuk, Trybala, and Starov (2014) N. Kovalchuk, A. Trybala, and V. Starov, Curr. Opin. Colloid Interface Sci. 19, 336 (2014).
  * Erbil and Meric (1997) H. Y. Erbil and R. A. Meric, J. Phys. Chem. B 101, 6867 (1997).
  * McHale, Newton, and Shirtcliffe (2009) G. McHale, M. I. Newton, and N. J. Shirtcliffe, J. Phys.: Condens. Matter 21, 464122 (2009).
  * Courbin  _et al._ (2007) L. Courbin, E. Denieul, E. Dressaire, M. Roper, A. Ajdari, and H. A. Stone, Nat Mater 6, 661 (2007).
  * Savva and Kalliadasis (2012) N. Savva and S. Kalliadasis, J. Eng. Math. 73, 3 (2012).
  * Modaressi and Garnier (2002) H. Modaressi and G. Garnier, Langmuir 18, 642 (2002).
  * Bliznyuk  _et al._ (2009) O. Bliznyuk, E. Vereshchagina, E. S. Kooij, and B. Poelsema, Phys Rev E 79, 041601 (2009).
  * Herminghaus, Brinkmann, and Seemann (2008) S. Herminghaus, M. Brinkmann, and R. Seemann, Annu. Rev. Mater. Res. 38, 101 (2008).
  * Tanner (1979) L. Tanner, J. Phys. D: Appl. Phys. 12, 1473 (1979).
  * Frank and Perre (2012) X. Frank and P. Perre, Phys. Fluids 24, 042101 (2012).
  * Acton, Huppert, and Worster (2001) J. M. Acton, H. E. Huppert, and M. G. Worster, J. Fluid Mech. 440, 359 (2001).
  * Starov (2004) V. M. Starov, Adv. Colloid Interface Sci. 111, 3 (2004).
  * He and Hadjiconstantinou (2003) G. He and N. Hadjiconstantinou, J. Fluid Mech. 497, 123 (2003).
  * De Ruijter, Blake, and De Coninck (1999) M. J. De Ruijter, T. Blake, and J. De Coninck, Langmuir 15, 7836 (1999).
  * Song, Li, and Liu (2013) F. Song, B. Li, and C. Liu, Langmuir 29, 4266 (2013).
  * Raiskinmäki  _et al._ (2000) P. Raiskinmäki, A. Koponen, J. Merikoski, and J. Timonen, Comput. Mater. Sci. 18, 7 (2000).
  * Attar and Körner (2009) E. Attar and C. Körner, J Colloid Interface Sci 335, 84 (2009).
  * Dupuis and Yeomans (2004) A. Dupuis and J. M. Yeomans, Future Generation Computer Systems 20, 993 (2004).
  * Ehrhard and Davis (1991) P. Ehrhard and S. H. Davis, J. Fluid Mech. 229, 365 (1991).
  * Espín and Kumar (2015) L. Espín and S. Kumar, J. Fluid Mech. 784, 465 (2015).
  * Zadražil, Stepanek, and Matar (2006) A. Zadražil, F. Stepanek, and O. Matar, J. Fluid Mech. 562, 1 (2006).
  * Greenspan (1978) H. P. Greenspan, J. Fluid Mech. 84, 125 (1978).
  * Cox (1986) R. Cox, J. Fluid Mech. 168, 169 (1986).
  * Haley and Miksis (1991) P. J. Haley and M. J. Miksis, J Fluid Mech 223, 57 (1991).
  * Bostwick (2013) J. Bostwick, J. Fluid Mech. 725, 566 (2013).
  * Smith (1995) M. K. Smith, J. Fluid Mech. 294, 209 (1995).
  * Gomba and Homsy (2010) J. M. Gomba and G. M. Homsy, J. Fluid Mech. 647, 125 (2010).
  * Hardy (1919) W. B. Hardy, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 38, 49 (1919).
  * Bascom, Cottington, and Singleterry (1964) W. D. Bascom, R. L. Cottington, and C. R. Singleterry (ACS Publications, 1964).
  * Radigan  _et al._ (1974) W. Radigan, H. Ghiradella, H. Frisch, H. Schonhorn, and T. Kwei, J. Colloid Interface Sci. 49, 241 (1974).
  * Ghiradella, Radigan, and Frisch (1975) H. Ghiradella, W. Radigan, and H. Frisch, J. Colloid Interface Sci. 51, 522 (1975).
  * Schwartz (1998) L. W. Schwartz, Langmuir 14, 3440 (1998).
  * Alleborn and Raszillier (2004) N. Alleborn and H. Raszillier, Chem. Eng. Sci. 59, 2071 (2004).
  * Davis and Hocking (1999) S. Davis and L. Hocking, Phys Fluids 11, 48 (1999).
  * De Ruijter  _et al._ (1997) M. J. De Ruijter, J. De Coninck, T. Blake, A. Clarke, and A. Rankin, Langmuir 13, 7293 (1997).
  * Kumar and Deshpande (2006) S. M. Kumar and A. P. Deshpande, Colloids Surf., A 277, 157 (2006).
  * Wang  _et al._ (2016) Z. Wang, L. Espín, F. S. Bates, S. Kumar, and C. W. Macosko, Chem. Eng. Sci. 146, 104 (2016).
  * Schwartz, Roux, and Cooper-White (2005) L. Schwartz, D. Roux, and J. Cooper-White, Physica D 209, 236 (2005).
  * Glasner (2003) K. Glasner, Phys Fluids 15, 1837 (2003).
  * Popescu  _et al._ (2012) M. N. Popescu, G. Oshanin, S. Dietrich, and A. Cazabat, J. Phys.: Condens. Matter 24, 243102 (2012).
  * Lafuma and Quéré (2003) A. Lafuma and D. Quéré, Nat Mater 2, 457 (2003).
  * Chow (1998) T. Chow, J. Phys.: Condens. Matter 10, L445 (1998).
  * Savva, Kalliadasis, and Pavliotis (2010) N. Savva, S. Kalliadasis, and G. A. Pavliotis, Phys Rev Lett 104, 084501 (2010).
  * Herminghaus (2012) S. Herminghaus, Eur. Phys. J. E 35, 1 (2012).


# Markovianity and non-Markovianity of Particle Bath with Dirac Dispersion
Relation

Takano Taira, Naomichi Hatano, Akinori Nishino

###### Abstract

The decay rate of quantum particles in open quantum systems has traditionally
been known as exponential, based on empirical predictions from experiments and
theoretical predictions from the Markovian dynamics of the corresponding
quantum states. However, both theoretical predictions and experimental
observations suggest deviations from this exponential decay, particularly in
the short and long time regimes. In this study, we explore the spontaneous
emission of a single Dirac particle within an environment characterized by an
energy spectrum with a gap $m$ and an energy cutoff $L$. Our results reveal
that high-energy structures, such as the spectral cutoff $L$, play a critical
role in driving the short-time non-exponential decay. In contrast, the long-
time decay is predominantly influenced by low-energy structures, such as the
Dirac gap $m$. Surprisingly, we find that in the limits where the energy
cutoff $L$ is infinite and the energy gap $m$ is zero, the decay dynamics of
massless Dirac particles exhibit Markovian behavior without the need for
conventional approximations like the Born-Markov approximation. This work
underscores the complex interplay between particle energy properties and decay
dynamics, providing new insights into quantum decay processes.

##  1 Introduction

The history of research on non-exponential decay in quantum systems dates back
to the work in 1958 by Khalfin [9] in the context of nuclear decay, which
predates the development of quantum master equations for open quantum systems,
such as the Nakajima-Zwanzig [15, 22] and
Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) [6, 11] equations, which were
developed in the 1960s and 1970s, respectively. Khalfin demonstrated a power-
law decay in the long-time regime, which was experimentally confirmed as late
as 2006 [20]. In 1977, Chiu, Sudarshan, and Misra [2] demonstrated a non-
exponential decay in the short-time regime, which led to the prediction of the
quantum Zeno effect [14], in which continuous observation of an excited state
prevents the state from decaying. This was experimentally observed in 1990
[8]. In these findings, Khalfin [9], Chiu, Misra, and Sudarshan [2] utilized a
theorem based on the Fourier analysis called the Paley-Wiener theorem [17],
which implies that if the spectrum of the Hamiltonian occupies the entire real
axis, the decay is purely exponential, and vice versa. This implies that when
the spectrum of the Hamiltonian is bounded from below, the decay profile is
non-exponential.

In a physical setup, it is natural to assume that the spectrum of the
Hamiltonian is bounded from below; hence, according to the prediction by
Khalfin, Chiu, Sudarshan, and Misra, the decay should be non-exponential in
general. However, the exponential decay is widely accepted as a standard decay
profile of quantum systems, as early as in the first empirical confirmation of
an exponential decay by Rutherford and Soddy in 1902 [21].

The exponential decay also follows from the GKSL equation. To see this,
consider the dynamical map of the system density
$\rho_{S}(t)=\exp(\mathcal{\mathcal{L}}t)\rho_{S}(0)$, where $\rho_{S}(t)$ is
the system density matrix at time $t$ and $\mathcal{L}$ is the Liouvillian
superoperator of the GKSL form. The decaying profile, or the survival
probability of an excited state $\ket{1}$ is given by
$\bra{1}\rho_{S}(t)\ket{1}$. If $\mathcal{L}$ is diagonalizable, the survival
probability decays exponentially with a possible periodic oscillation. When
$\mathcal{L}$ is non-diagonalizable due to an exceptional point [13], the
survival probability initially shows a power-law increase but is later
suppressed by an exponential decay. In short, the GKSL equation shows an
exponential decay. Since GKSL (or GKSL-type) equation follows from the
assumption of semi-group property of the dynamical map, often taken as the
definition of the quantum Markovianity [19], one can deduce that the Markovian
dynamics of the system implies an exponential decay. The negation of such a
statement is that the non-exponential decay implies non-Markovian dynamics of
the system. Therefore, in a physical setup, in which the spectrum of the full
Hamiltonian is bounded from below, the system’s dynamics must be non-
Markovian, and only with Markov and Born approximations can the Markovian
dynamics be achieved.

Short-time decay profile of survival prob. $m=0$ $m\not=0$ $L<\infty$
Quadratic Quadratic $L\rightarrow\infty$ Exp. Exp.

Long-time decay profile of survival prob. $m=0$ $m\not=0$ $L<\infty$ Power law
decay to bound state Power law decay to bound state $L\rightarrow\infty$ Exp.
decay to zero Power law decay to bound state

Table 1: Two tables showing the decay profiles of the survival probability
$P(t)$ for short and long time regimes. The picture shows a schematic of our
model (1) with coupling $g$, gap $m$, and momentum cut off $L$ shown.

On the other hand, it has been demonstrated analytically in electron transport
in a single-electron case [7] and an interacting case [16] that the excited
state decays exponentially without any approximation. This phenomenon can be
attributed to the model being understood as a massless Dirac particle, whose
spectrum features a Dirac cone and is unbounded. Thus, the question arises:
What happens to the decay profile when a Dirac gap opens? This is the central
question of this paper.

We will consider a simple two-level system attached to a Dirac-particle bath,
schematically represented in Table 1. Our result is also summarised in Table
1, where high-energy structures, such as the cutoff $L$, primarily contribute
to the short-time regime producing a quadratic decay, whereas low-energy
structures, such as the Dirac gap $m$, mainly contribute to the long-time
regime producing a power-law decay and localized bound states, as seen in Ref.
[10, 3]. Surprisingly, the opening of the Dirac gap, which has no lower or
upper bound, still induces deviations from pure exponential decay. We note
that the correspondence between the gap structure of the spectrum and the
survival probability was also studied in the context of quantum optics [10]
and with the tight-binding model [3] in long-time regimes. In the present
work, we have extended their work to the Dirac gap $m$, and further analyzed
the interplay with the cutoff $L$ in both short- and ling-time regimes.

Our model is mathematically similar to quantum optical models. Two key
differences are that we consider a Dirac particle instead of photons and that
we focus on controlling the Dirac gap to observe various non-Markovian
dynamics; in quantum optics, on the other hand, the Dirac cone is typically
fixed [18, 5].

Our model has a possibility for applicaitons in photonic crystals, which can
be engineered to exhibit a linear frequency dispersion [18, 5]. In this case,
the spectrum of the Hamiltonian is unbounded, with negative energy states
occupied by antiparticles, or holes. Therefore, our results in Table 1 suggest
that such systems might exhibit Markovian dynamics without any approximations.
Moreover, when the mass of the Dirac particle is finite, a gap opens in the
Hamiltonian’s spectrum, potentially exhibiting non-Markovian dynamics in the
long-time regime. These insights suggest that controlling the spectrum’s gap
to regulate the Markovianity of system dynamics could benefit both quantum
optics and electron systems.

In Section 2, we introduce our toy model and derive the master equation. In
Section 3 we analytically investigate the master equation in the cases of
$m>0,L\rightarrow\infty$ and $m=0,L<\infty$. In Section 4, we numerically
investigate the most general case $m>0,L<\infty$.

##  2 Survivial probability of two-level system with half-filled Dirac
environment

In this section, we will introduce the model and obtain a non-Markovian master
equation of the system, starting from the Schrödinger equation of the full
system.

###  2.1 Model

We introduce a simple toy model, designed to capture the essence of the
interplay between the Dirac particle’s mass and the resulting dynamical
behavior:

$\displaystyle
H_{\text{full}}=\omega_{0}\sigma_{+}\sigma_{-}+H_{D}+\left(\sigma_{+}\otimes
B+\sigma_{-}\otimes B^{\dagger}\right).$ (1)

The first term of the full Hamiltonian (1) is the two-level Hamiltonian, which
we refer to as the two-level system with the states
$\\{\ket{0}_{S},\ket{1}_{S}\\}$. We take the level difference $\omega_{0}$ to
be real and constant. The operators $\sigma_{+}:=(\sigma_{x}-i\sigma_{y})/2$
and $\sigma_{-}:=(\sigma_{x}+i\sigma_{y})/2$ are the creation and annihilation
operators of the two-level system, which satisfies the fermionic anti-
commutation relations:

$\displaystyle\\{\sigma_{-},\sigma_{+}\\}=\\{\sigma_{-},\sigma_{+}\\}=1,\quad\\{\sigma_{-},\sigma_{-}\\}=\\{\sigma_{+},\sigma_{+}\\}=0.$
(2)

The second term of Eq. (1) is the half-filled Dirac Hamiltonian, which we
refer to as the enviroment:

$\displaystyle:H_{D}:=\int_{-L}^{L}dp~{}\omega_{p}\left(b^{\dagger}_{p}b_{p}-c_{p}c^{\dagger}_{p}\right)\text{
with }\quad\ket{\text{vac}}_{B}:=\prod_{p}c^{\dagger}_{p}\ket{0},$ (3)

where $L>0$ is the cutoff on the momentum space. We note that the normal
ordering is taken with respect to the vacuum state $\ket{\text{vac}}_{B}$. The
operators $b^{\dagger}_{p}$ and $c_{p}^{\dagger}$ are creation operators of
electrons with the same dispersion relation $\omega_{p}$, satisfying the
commutation relations:

$\displaystyle\\{b_{p},b^{\dagger}_{k}\\}=\\{c_{p},c^{\dagger}_{k}\\}=\delta(p-k),$
(4)
$\displaystyle\\{b_{p},b_{k}\\}=\\{c_{p},c_{k}\\}=\\{b_{p},c_{k}\\}=\\{b_{p},c_{k}^{\dagger}\\}=0,$
(5)

where curly brackets indicate the anti-commutator. The dispersion relation for
massive and massless particles is given by

massive $\displaystyle:$ $\displaystyle\omega_{p}=\sqrt{p^{2}+m^{2}},$ (6)
massless $\displaystyle:$ $\displaystyle\omega_{p}=|p|.$ (7)

Finally, the last two terms of Eq. (1) represent the interaction between the
two-level system and the environment. The third term $\sigma_{+}\otimes B$,
where $B:=\int_{-L}^{L}dp~{}g_{p}B_{p}$, $B_{p}:=b_{p}+c_{p}^{\dagger}$,
creates the particle in the two-level system and at the same time either
destroys a particle by $b_{p}$ or creates a particle by $c_{p}^{\dagger}$ in
the environment, preserving the overall particle number. The coupling
parameter $g_{p}$ captures the details of the interaction. Throughout this
paper, we consider the case of a simple uniform coupling $g_{p}=g$.

###  2.2 Master equation

Next, let us find the master equation associated with our model. For brevity,
let us consider the following abbreviation of the composite states:

$\displaystyle\ket{0\text{vac}}$ $\displaystyle:=$
$\displaystyle\ket{0}_{S}\otimes\ket{\text{vac}}_{B},$ (8)
$\displaystyle\ket{1\text{vac}}$ $\displaystyle:=$
$\displaystyle\ket{1}_{S}\otimes\ket{\text{vac}}_{B}.$ (9)

Note that $\ket{1\text{vac}}=\sigma^{+}\ket{0\text{vac}}$. Following the same
procedure as in Refs. [4, 1], we let the initial state have a single particle
in the excited state of the two-level system:

$\displaystyle\ket{\phi(0)}:=\phi(0)\ket{1\text{vac}}.$ (10)

Given that the full system (1) is assumed to be closed, the particle number is
conserved. This can be checked by calculating the commutation relation between
the Hamiltonian (1) and the number operator. Therefore, the state at an
arbitary time $t$ must be given by

$\displaystyle\ket{\phi(t)}:=\phi(t)\ket{1\text{vac}}+\int_{-L}^{L}dp~{}\phi_{p}(t)b_{p}^{\dagger}\ket{0\text{vac}}+\int_{-L}^{L}dp~{}\psi_{p}(t)c_{p}\ket{0\text{vac}}.$
(11)

We require the coefficients to satisfy the normalization condition:

$\displaystyle\left|\phi(t)\right|^{2}+\int_{-L}^{L}dp~{}\left|\phi_{p}(t)\right|^{2}+\int_{-L}^{L}dp~{}\left|\psi_{p}(t)\right|^{2}=1.$
(12)

Inserting the Ansatz (11) to the Schrödinger equation of the full system
$i\hbar\partial_{t}\ket{\psi}=H_{\text{full}}\ket{\psi}$, we find the
following coupled first-order equations of the coefficients $\phi,\phi_{p}$
and $\psi_{p}$:

$\displaystyle i\hbar\frac{d\phi(t)}{dt}$
$\displaystyle=\omega_{0}\phi(t)+g\int_{-L}^{L}dp~{}(\phi_{p}(t)+\psi_{p}(t)),$
(13) $\displaystyle i\hbar\frac{d\phi_{p}(t)}{dt}$
$\displaystyle=\omega_{p}\phi_{p}(t)+g\phi(t),$ (14) $\displaystyle
i\hbar\frac{d\psi_{p}(t)}{dt}$
$\displaystyle=-\omega_{p}\psi_{p}(t)+g\phi(t).$ (15)

In order to eliminate the first term on the right-hand side of each equation,
we introduce the new functions $\Phi(t):=e^{i\omega_{0}t/\hbar}\phi(t)$,
$\Phi_{p}(t):=e^{i\omega_{p}t/\hbar}\phi_{p}(t)$ and
$\Psi_{p}(t):=e^{-i\omega_{p}t/\hbar}\psi_{p}(t)$. The coupled first-order
equations now take the following form:

$\displaystyle i\hbar\frac{d\Phi(t)}{dt}$
$\displaystyle=g\int_{-L}^{L}dp~{}e^{\frac{i}{\hbar}\omega_{0}t}\left(\Phi_{p}(t)e^{-\frac{i}{\hbar}\omega_{p}t}+\Psi_{p}(t)e^{\frac{i}{\hbar}\omega_{p}t}\right),$
(16) $\displaystyle i\hbar\frac{d\Phi_{p}(t)}{dt}$
$\displaystyle=ge^{-\frac{i}{\hbar}(\omega_{0}-\omega_{p})t}\Phi(t),$ (17)
$\displaystyle i\hbar\frac{d\Psi_{p}(t)}{dt}$
$\displaystyle=ge^{-\frac{i}{\hbar}(\omega_{0}+\omega_{p})t}\Phi(t).$ (18)

The last two equations can be solved, yielding the following equations:

$\displaystyle\Phi_{p}(t)$
$\displaystyle=\frac{1}{i\hbar}g\int_{0}^{t}ds~{}e^{-\frac{i}{\hbar}(\omega_{0}-\omega_{p})s}\Phi(s)+\Phi_{p}(0),$
(19) $\displaystyle\Psi_{p}(t)$
$\displaystyle=\frac{1}{i\hbar}g\int_{0}^{t}ds~{}e^{-\frac{i}{\hbar}(\omega_{0}+\omega_{p})s}\Phi(s)+\Psi_{p}(0).$
(20)

Inserting these to Eq. (16), the explicit form of the coefficient $\Phi(t)$ is
found by solving the following non-Markovian master equation:

$\displaystyle\frac{d\Phi(t)}{dt}$ $\displaystyle=$
$\displaystyle-\frac{1}{\hbar^{2}}\int^{t}_{0}dse^{\frac{i}{\hbar}\omega_{0}(t-s)}\Phi(s){g^{2}\int_{-L}^{L}dp~{}\cos\left[\omega_{p}\frac{t-s}{\hbar}\right]}$
(21)
$\displaystyle-\frac{i}{\hbar}g^{2}\int_{-L}^{L}dp~{}e^{\frac{i}{\hbar}\omega_{0}t}\left(e^{-\frac{i}{\hbar}\omega_{p}t}\Phi_{p}(0)+e^{\frac{i}{\hbar}\omega_{p}t}\Psi_{p}(0)\right).$

Notice that with our initial condition (10), the second term of the above
equation vanishes. We will refer to the $p$ integral part of the first term as
the memory kernel.

The norm squared of the solution to the above master equation is the survival
probability of the particle at the state $\ket{1\text{vac}}$. To see this, we
will follow the procedure in Sections 10.1.1 and 10.1.2 of Ref. [1]. We first
define the system density matrix as

$\displaystyle\rho_{S}(t):=\text{Tr}_{B}\left(\ket{\phi(t)}\bra{\phi(t)}\right),$
(22)

where the trace over $B$ is defined as

$\displaystyle\text{Tr}_{B}\left({\rho}\right):={{}_{B}\bra{\text{vac}}}{\rho}\ket{\text{vac}}_{B}+\sum_{p}{{}_{B}\bra{\text{vac}}b_{p}}{\rho}b_{p}^{\dagger}\ket{\text{vac}}_{B}+\sum_{p}{{}_{B}\bra{\text{vac}}c_{p}^{\dagger}}{\rho}c_{p}\ket{\text{vac}}_{B}.$
(23)

Using the normalization condition (12), we arrive at the system density matrix
of the form

$\displaystyle\rho_{S}(t)$ $\displaystyle=$
$\displaystyle\left(1-\left|\phi(t)\right|^{2}\right)\ket{0}_{SS}\bra{0}+\left|\phi(t)\right|^{2}\ket{1}_{SS}\bra{1}$
(29) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\ket{0}_{S}&\ket{1}_{S}\\\
\end{array}\right)\left(\begin{array}[]{cc}1-\left|\phi(t)\right|^{2}&0\\\
0&\left|\phi(t)\right|^{2}\end{array}\right)\left(\begin{array}[]{c}{{}_{S}\bra{0}}\\\
{{}_{S}\bra{1}}\end{array}\right).$

The survival probability of the qubit state $\ket{1}_{S}$ is therefore given
by

$\displaystyle
P(t):=\left|\bra{1}\rho_{S}(t)\ket{1}\right|=\left|\phi(t)\right|^{2}=\left|\Phi(t)\right|^{2}.$
(30)

The last equality follows from the definition that $\phi$ and $\Phi$ are
different only up to the phase. This quantity’s physical meaning is the
particle’s probability at the state $\ket{1}_{S}$ to remain or return to the
original state $\ket{1}_{S}$ after time $t$.

To solve the master equation (21), we will utilize the Laplace transform.
First, for brevity, let us introduce the memory kernel $\mathcal{K}(t)$ to
rewrite the master equation (21):

$\displaystyle\frac{d\Phi(t)}{dt}=-\frac{g^{2}}{\hbar^{2}}\int_{0}^{t}ds~{}\Phi(s)\mathcal{K}(t-s),\quad\mathcal{K}(\tau):=e^{\frac{i}{\hbar}\omega_{0}\tau}\int_{-L}^{L}dp~{}\cos\left(\omega_{p}\frac{\tau}{\hbar}\right).$
(31)

After applying the Laplace transform detailed in Appendix C, we find the
formal solution to the master equation (21):

$\displaystyle\Phi(t)=e^{\frac{i}{\hbar}\omega_{0}t}\lim_{R\rightarrow\infty}\int_{-iR+\sigma}^{iR+\sigma}\frac{dz}{2\pi
i}e^{z\frac{t}{\hbar}}\left(z+i\omega_{0}+2g^{2}\frac{z\text{Arctan}\left(\frac{L}{\sqrt{m^{2}+z^{2}}}\right)}{\sqrt{m^{2}+z^{2}}}\right)^{-1}\\!\\!\\!\\!\\!\\!\Phi(0).$
(32)

##  3 Spectural structure and Markovianity

The integral in the solution (32) can be evaluated numerically. It shows a
$t^{-2}$ decay in a short time, an exponential decay in mid-time, and a power
law decay with oscillation in a long time. In this section, we investigate the
relation between these non-exponential behaviors and the spectrum structure.
We will find that the low-energy structure, such as the energy gap,
corresponds to the long-time behavior, whereas the high-energy structure, such
as the cutoff, corresponds to the short-time behavior; see Table 1 and for the
long-time behavior in Table 1, where, lightly shaded regions indicate an
exponential decay, while darkly shaded regions indicate a non-exponential
decay.

###  3.1 Low energy structure and Markovianity: $m>0,L\rightarrow\infty$ limit

First, consider the limit $L\rightarrow\infty$, which removes the high-energy
cutoff and leaves the low-energy energy gap $2m$. In this case, the Arctan
function in the integrand of Eq. (32) becomes $\pi/2$, and hence the solution
(32) is reduced to the following relatively simple form:

$\displaystyle\Phi(t)$ $\displaystyle=$ $\displaystyle
e^{\frac{i}{\hbar}\omega_{0}t}\lim_{R\rightarrow\infty}\int_{-iR+\sigma}^{iR+\sigma}\frac{dz}{2\pi
i}e^{z\frac{t}{\hbar}}F(z)\Phi(0),$ (33)

where we have defined the complex function $F(z)$ as

$\displaystyle F(z):=\left(z+i\omega_{0}+g^{2}\frac{\pi
z}{\sqrt{m^{2}+z^{2}}}\right)^{-1}.$ (34)

Notice that the function $F(z)$ depends on the square root of the complex
variable $z$, which can only be smoothly defined on the Riemann surface. Let
us denote the two Riemann sheets as $\mathcal{R}_{\pm}$, where labeling
indicates the sign $\sqrt{-1}=\pm i$. We draw a branch cut from the point
$z=\pm im$ going off to negative infinity as shown in Fig. 1(a). Note that the
infinity is not the branch points of $\sqrt{m^{2}+z^{2}}$.

Figure 1: Panel (a) shows the contour, branch cuts, and poles of the Bromwich
integral (33), which is equivalent to the contour shown in panel (b).

Typically, the above integral is computed using the Cauchy residue theorem,
simplifying the integral to the sum of the residues of the integrand. However,
due to the branch cut in the function $F(z)$, not all poles are located on the
first Riemann sheet, preventing a straightforward summation of all pole
contributions.

In order to find the correct locations of the poles, we reexpress the function
as follows:

$\displaystyle F(z)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{m^{2}+z^{2}}}{\left(z+i\omega_{0}\right)\sqrt{m^{2}+z^{2}}+g^{2}\pi
z}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{m^{2}+z^{2}}}{\left(z+i\omega_{0}\right)\sqrt{m^{2}+z^{2}}+g^{2}\pi
z}\frac{\left(z+i\omega_{0}\right)\sqrt{m^{2}+z^{2}}-g^{2}\pi
z}{\left(z+i\omega_{0}\right)\sqrt{m^{2}+z^{2}}-g^{2}\pi z}$ $\displaystyle=$
$\displaystyle\frac{\left(z+i\omega_{0}\right)(m^{2}+z^{2})-g^{2}\pi
z\sqrt{m^{2}+z^{2}}}{\mathcal{P}(z)},$ (36)

where $\mathcal{P}(z)=(z+i\omega_{0})^{2}(m^{2}+z^{2})-g^{4}\pi^{2}z^{2}$
represents the fourth-order polynomial in the denominator of the final
expression. The zeros of $\mathcal{P}(z)$ determine the four poles of the
integrand. Analyzing the discriminant of the polynomial, which can be
evaluated explicitly, reveals the presence of two purely imaginary solutions,
denoted by $z_{0},z_{1}\in i\mathbb{R}$, and a pair of complex solutions,
$z_{\pm}$, where $z_{+}=-(z_{-})^{*}$. Imposing the constraint that parameters
$\omega_{0},m,g$ are all positive, one can show explicitly that the imaginary
parts of $z_{0},z_{1}$ are negative and positive, respectively. Notice that
$z_{\pm}$ may also be purely imaginary under certain conditions. The regime
allowing for purely imaginary poles $z_{\pm}$ can be deduced by examining the
discriminant of $\mathcal{P}=0$, where the phase diagram, illustrated in Fig.
2, depicts a transition from four imaginary poles to two imaginary and two
complex poles. This transition point is known as the exceptional point.

Going back to the expression of the integrand (3.1), one can show by using the
fact that on the principle Riemann sheet $\mathcal{R}_{+}$, the square root is
defined as $\sqrt{-1}=i$, only $z_{0}$ is the pole of the function $F(z)$.
Other zeros of $\mathcal{P}(z)$ are poles of $F(z)$ only on the second Riemann
sheet $\mathcal{R}_{-}$, where $\sqrt{-1}=-i$:

$\begin{array}[]{cc}z_{0}\in-i\mathbb{R}_{\geq 0}&\text{ on
}\mathcal{R}_{+},\\\ z_{1}\in i\mathbb{R}_{\geq 0},\quad
z_{+}=-(z_{-})^{*}&\text{ on }\mathcal{R}_{-}.\end{array}$ (37)

This would lead to the numerator of Eq. (36) vanishing for zeros
$z_{1},z_{\pm}$. Note that contrary to the conventional Markovian dynamics,
where the exceptional point of the Lindbladian operator shows a sudden change
in the dynamical behavior of the state, in our case, the exceptional point in
the parameter space shown in Fig. 2 does not affect to the system’s dynamics.

Figure 2: Phase diagram of solutions to the fourth-order polynomial
$\mathcal{P}(z)=(z+i\omega_{0})^{2}(m^{2}+z^{2})-g^{4}\pi^{2}z^{2}$, where we
put $\hbar=\omega_{0}=1$.

We are now in a position to evaluate the Bromwich integral in Eq. (33), where
its contour and locations of the poles and branch cuts are shown in Fig. 1(a).
The Bromwich contour in Eq. (33), shown in Fig. 1(a) is equivalent to the
contour in Fig. 1(b). This can be shown by using the Cauchy theorem; see
details in Appendix B. Therefore the Bromwich integral in Eq. (33) is
decomposed to the contributions from the pole $z_{0}$ on the first Riemann
sheet and the contour around the branch cut:

$\displaystyle\Phi(t)=e^{\frac{i}{\hbar}\omega_{0}t}\left[\oint_{z_{0}}\dots+\int_{\text{BC}}\dots\right]\Phi(0).$
(38)

Applying the Cauchy residue theorem at the pole $z_{0}$, we obtain

$\displaystyle\oint_{z_{0}}\dots=\text{Res}(F(z)e^{zt},z_{0})=\frac{\left(z_{0}+i\omega_{0}\right)(m^{2}+z^{2}_{0})-g^{2}\pi
z_{0}\sqrt{m^{2}+z^{2}_{0}}}{(z_{0}-z_{1})(z_{0}-z_{+})(z_{0}-z_{-})}e^{z_{0}t}$
(39)

Given that $z_{0}$ is purely imaginary, this value does not contribute to the
system’s decay but to a bound state.

Note that if we set $\omega_{0}=0$ and evaluate the residue, we find the
following simple expression:

$\displaystyle\oint_{z_{0}}\dots=\text{Res}(F(z)e^{zt},z_{0})=\frac{m}{m+\pi
g^{2}}.$ (40)

Note that this expression cannot be obtained by simply taking the limit
$\omega_{0}\rightarrow 0$ in Eq. (39), as one must take the limit there is a
limit $\omega_{0}\rightarrow 0$, before calculating the residue.

Next, let us consider the branch cut contribution. Let $F_{\pm}(z)$ denote the
integrand of Eq. (33) on the Riemann sheets $\mathcal{R}_{\pm}$, respectively.
We explicitly have

$\displaystyle
F_{\pm}(z)=\frac{\left(z+i\omega_{0}\right)(m^{2}+z^{2})}{\mathcal{P}(z)}\pm\frac{-g^{2}\pi
z\sqrt{m^{2}+z^{2}}}{\mathcal{P}(z)},$ (41)

where the square-root function in the above expression is assumed to be the
principal branch with $\sqrt{-1}=i$.

The upper part of the branch-cut contribution in Fig. 1(b) is given by

$\displaystyle\int_{-R+im}^{im}\frac{dz}{2\pi
i}~{}F_{+}(z)e^{z\frac{t}{\hbar}}+\int^{-R+im}_{im}\frac{dz}{2\pi
i}~{}F_{-}(z)e^{z\frac{t}{\hbar}}$ (42) $\displaystyle=$
$\displaystyle\int_{-R+im}^{im}\frac{dz}{2\pi
i}~{}\left(F_{+}(z)-F_{-}(z)\right)e^{z\frac{t}{\hbar}}$ $\displaystyle=$
$\displaystyle\int_{-R+im}^{im}\frac{dz}{2\pi i}~{}A(z)e^{z\frac{t}{\hbar}}$
$\displaystyle=$ $\displaystyle\int_{-R}^{0}\frac{dw}{2\pi
i}~{}A\left(w+im\right)e^{w\frac{t}{\hbar}}e^{im\frac{t}{\hbar}},$ (43)

and $A(z):=(-2g^{2}\pi z\sqrt{m^{2}+z^{2}})/\mathcal{P}(z)$, where we have
performed a coordinate transformation $z=w+im$. The lower part of the branch
cut contribution can be calculated in a similar way, and thus we find the full
branch-cut contribution in the form

$\displaystyle\int_{\text{BC}}\dots=\lim_{R\rightarrow\infty}\int_{-R}^{0}\frac{dw}{2\pi
i}~{}\left(A\left(w+im\right)e^{im\frac{t}{\hbar}}-A\left(w-im\right)e^{-im\frac{t}{\hbar}}\right)e^{w\frac{t}{\hbar}}.$
(44)

For the rest of this paper, we will refer to Eq. (39) and the integral form of
Eq. (44) inserted into Eq. (38) as the ”exact solution.”

Given explicit expression of the integral, we turn to the large-time
asymptotic analysis outlined in Ref. [3]. First, we change the complex varible
from $\omega$ to $\tau$, related via $wt\equiv\tau$. The integrand can be
expanded as a series of $1/t$. The lowest-order term is $1/\sqrt{t}$, which
after evaluation gives

$\displaystyle\int_{\text{BC}}\dots=\frac{1}{g^{2}\sqrt{m}}\left(\cos\left(\frac{mt}{\hbar}\right)+\sin\left(\frac{mt}{\hbar}\right)\right)\left(\frac{\hbar}{\pi
t}\right)^{3/2}+\mathcal{O}\left(\frac{1}{t^{5/2}}\right).$ (45)

Therefore, we find the asymptotic contribution of a branch cut to be the
power-law decay.

In Fig. 3, we present our findings. Figure 3 illustrates the survival
probability $P(t)=\left|\Phi(t)\right|^{2}$ over the entire time range. We
note that $P(t)$ decays to a constant value over time, which is the residue of
the pole $z_{0}$:

$\displaystyle\left|\text{Res}(F(z)e^{zt},z_{0})\Phi(0)\right|^{2}.$ (46)

Given that the pole $z_{0}$ is always pure imaginary, the magnitude of this
residue is independent of time, allowing the survival probability to maintain
a finite value indefinitely. The constant value predicted by the residue is
depicted as a horizontal red dashed line in Fig. 3(a).

The short-time dynamics are depicted in Fig. 3(b), where we used a log-log
plot to display $1-P(t)$ for the exact solution (38). The slope of the
straight line on this plot represents the power governing the decay. From Fig.
3(b), it is evident that the exact solution exhibits a linear decay at short
times. This observation suggests that the quantum Zeno effect (see Appendix
A), which requires a quartic short-time decay, is absent when
$L\rightarrow\infty$. Consequently, we infer that the quantum Zeno effect
arises from the spectral bounds of the system rather than the gap structure.
We will see in the next sub section 3.2 that a finite value of $L$ introduces
a new pair of poles, which gives rise to the quartic decay.

Finally, to see the decay in the long-time regime, we subtract the bound-state
contribution from the exact solution (38) and plot only the branch cut
contribution in Fig. 3(c). The dahsed lines in Fig. 3(c) indicates the
approximations upto order $t^{-3/2}$ and upto $t^{-5/2}$ given by Eq. (45). We
observe a good agreement with the exact solution beyond $t=10$, and from Fig.
3(c), we observe that the approximation gives a good prediction beyond $t=10$.

Figure 3: All three panels are showing the plot for
$\hbar=1,\omega_{0}=0.1,m=1,g=0.6$ in different time regimes. (a) the plot of
the survival probability $P(t)$ given by Eq. (38). The horizontal red dashed
line indicates the absolute value squared of the pole contribution (39) (b)
the log-log plot of $1-P(t)$ for the exact solution (38), plotted along $t$
and $t^{2}$. (c) the log plot of the numerical solution and approximations up
to $t^{-3/2}$and upto $t^{-5/2}$ in Eq. (45).

###  3.2 High energy structure and Markovianity: $m=0,L<\infty$ case

To investigate the high-energy contribution for $L<\infty$, let us remove the
low-energy contribution by taking $m=0$. The formal solution is simplified to

$\displaystyle\Phi(t)$ $\displaystyle=$ $\displaystyle
e^{\frac{i}{\hbar}\omega_{0}t}\lim_{R\rightarrow\infty}\int_{-iR+\sigma}^{iR+\sigma}\frac{dz}{2\pi
i}e^{z\frac{t}{\hbar}}\left(z+i\omega_{0}+2g^{2}\text{Arctan}\left(\frac{L}{z}\right)\right)^{-1}\Phi(0)$
(47) $\displaystyle=$ $\displaystyle
e^{\frac{i}{\hbar}\omega_{0}t}\lim_{R\rightarrow\infty}\int_{-iR+\sigma}^{iR+\sigma}\frac{dz}{2\pi
i}e^{\frac{1}{\hbar}g^{2}zt}\left(z+ia+2\text{Arctan}\left(\frac{b}{z}\right)\right)^{-1}\Phi(0),$

where in the last line, we have rescaled the integration variable from $z$ to
$z/g^{2}$ and introduced new parameters $a\equiv\omega_{0}/g^{2}$ and $b\equiv
L/g^{2}$.

The integrand now has a different type of Riemann surface, where branch points
are located at $z=\pm ib=\pm iL/g^{2}$ and $\infty$. To find the poles, we
must find the solutions to the following equation:

$\displaystyle
z+ia+2\text{Arctan}\left(\frac{b}{z}\right)=z+ia-i\text{Log}\left(\frac{z+ib}{z-ib}\right)=0.$
(48)

This transcendental equation cannot be solved directly. However, by defining
$w=-iz-b$, $r=-\exp(-a-b)$, and $\alpha=-2br$, we can rewrite it as

$\displaystyle we^{w}+rw=\alpha,\quad
r\in\left[-1,0\right),\quad\alpha\in\left(0,2b\right].$ (49)

This is a rewritten form of the generalized Euler-Lambert or $r$-Lambert
equation [12]. The solution to the above equation is denoted by the
$r$-Lambert function of the parameter $\alpha$: $W_{r}(\alpha)$. The branch
structure of this function has been rigorously analyzed in Ref.[12], revealing
that there are only two real solutions on the principal Riemann sheet, which
translate to purely imaginary solutions in the $z$ plane. While infinitely
many solutions exist for the above equation, only two real solutions can be
found on the principle Riemann sheet for $\alpha\in(0,2b)$.

From the above analysis, we find two pure imaginary poles $ix_{1}$ and
$-ix_{2}$, where $x_{1},x_{2}>b$, and branch cut of the integrand of Eq. (47).
Two examples of poles and branch points locations are shown in Fig.4(a) for
parameters $L=2,g=1,\omega_{0}=0$, and in Fig.4(b) for $L=2,g=1,\omega_{0}=8$.
Notice that locations of poles are almost in complex conjugate pairs in
Fig.4(a) and unbalanced in Fig.4(b).

Following the trick similar to the one in the previous section, we can rewrite
the Bromwich integral into the contributions from the poles and branch cut:

$\displaystyle\Phi(t)=e^{\frac{i}{\hbar}\omega_{0}t}\left[\oint_{ix_{1}}\dots+\oint_{-ix_{2}}\dots+\int_{\text{BC}}\dots\right]\Phi(0),$
(50)

where $ix_{1}$ and $-ix_{2}$ are two purely imaginary solutions to Eq. (48)
with $x_{1},x_{2}>b$.

Figure 4: Panel (a) and (b) show the plot of poles on the principle Riemann
surface for $L=2,g=1,\omega_{0}=0$ and $L=2,g=1,\omega_{0}=8$ respectively.

####  3.2.1 Pole contributions

Let us look into the contribution from the poles, starting with the upper
pole:

$\displaystyle\oint_{ix_{1}}\dots$ $\displaystyle=$
$\displaystyle\text{Res}\left[F(z)e^{zt},ix_{1}\right]=\lim_{z\rightarrow
ix_{1}}\frac{\left(z-ix_{1}\right)e^{zt}}{z+i\frac{\omega_{0}}{\hbar}+\frac{2g^{2}}{\hbar}\text{Arctan}\left(\frac{L}{z\hbar}\right)}.$
(51)

This limit is challenging to evaluate, as the denominator cannot be expressed
in polynomial form. To overcome this, we apply the l’Hôpital’s rule, which
equates the limit of a fraction to the limit of the fraction of derivatives,
thereby simplifying the calculation of the residue:

$\displaystyle\oint_{ix_{1}}\dots$ $\displaystyle=$
$\displaystyle\lim_{z\rightarrow
ix_{1}}\frac{\left(z-ix_{1}\right)e^{zt}}{z+i\frac{\omega_{0}}{\hbar}+\frac{2g^{2}}{\hbar}\text{Arctan}\left(\frac{L}{z\hbar}\right)}=\frac{\left(x_{1}-b\right)\left(b+x_{1}\right)e^{\frac{ig^{2}tx_{1}}{\hbar}}}{x_{1}^{2}-(b-2)b},$
(52)

where $b=L/g^{2}$.

Similarly, the contribution from the lower pole is determined by the
l’Hôpital’s rule again. We derive the combined contribution of both poles,
$\text{Pole}(t):=\oint_{ix_{1}}+\oint_{-ix_{2}}$ as

$\displaystyle\text{Pole}(t)=\frac{\left(x_{1}-b\right)\left(b+x_{1}\right)e^{i\frac{g^{2}x_{1}}{\hbar}t}}{x_{1}^{2}-(b-2)b}+\frac{\left(x_{2}-b\right)\left(b+x_{2}\right)e^{-i\frac{g^{2}x_{2}}{\hbar}t}}{x_{2}^{2}-(b-2)b}.$
(53)

Notice that when two poles $ix_{1}$ and $-ix_{2}$ are almost complex conjugate
to each other as seen in Fig.4(a), which means $x_{1}\sim x_{2}$, then the
above pole contribution is reduced to the periodic contribution:

$\displaystyle\text{Pole}(t)$ $\displaystyle\sim$
$\displaystyle\frac{2\left(x_{1}^{2}-b^{2}\right)\cos\left(\frac{g^{2}x_{1}}{\hbar}t\right)}{x_{1}^{2}-(b-2)b}.$
(54)

When one of the poles is close to the branch point but one is away from the
branch point as shown in the Fig.4(b), then the pole contribution gives bound-
state contribution

$\displaystyle\text{Pole}(t)$ $\displaystyle\sim$
$\displaystyle\frac{\left(x_{2}^{2}-b^{2}\right)e^{-i\frac{g^{2}x_{2}t}{\hbar}}}{x_{2}^{2}-(b-2)b}.$
(55)

Therefore, the pole contribution shows no decay process and only bound-state
or periodic contribution for different parameter values.

Figure 5: (a) Contour plot of $\text{exp}(-|x_{1}-x_{2}|)$, showing the
periodic regions where $x_{1}\sim x_{2}$ and bound region where $x_{1}$ and
$x_{2}$ are different. Two points indicate the parameter choices for the plot
of the survival probability in Fig.6. (b) Contour plot of
$-\log(|\text{Pole}(0)|^{2})$, with overlap of the contour lines of
$\text{exp}(-|x_{1}-x_{2}|)$.

The above analysis reveals that the periodicity and bound-state behavior of
the survival probability depend on the distances from each pole, $ix_{1}$ and
$-ix_{2}$, to the branch points $iL/g^{2}$ and $-iL/g^{2}$. Therefore, to
investigate the parameter regimes where periodic or bound-state behavior is
observed, we plot the contour of the function $e^{-|x_{1}-x_{2}|}$ over the
parameters $a=\omega_{0}/g^{2}$ and $b=L/g^{2}$ as shown in Fig.5(a). Notice
the smooth transition from the periodic to the bound-state contribution for
$(\omega_{0}/g^{2},L/g^{2})=(0,5)$ to $(\omega_{0}/g^{2},L/g^{2})=(10,5)$ with
corresponding survival probabilities shown in Fig.6(a) and Fig.6(b).
respectively.

Comparing two contour plots Fig.5(a) and Fig.5(b), one may notice that most of
the periodic region shown in Fig.5(a) coincide with the yellow to red region
of Fig.5(b). Only periodic region that is within $10^{0}$ to $10^{-5}$
amplitude in Fig.5(b) is for smaller values of $\omega_{0}/g^{2}$. As an
example, we have taken $(\omega_{0}/g^{2},L/g^{2})=(0,5)$, which indeed shows
dominating periodic contribution in Fig.6(a). By increasing $\omega_{0}/g^{2}$
to $(\omega_{0}/g^{2},L/g^{2})=(10,5)$, we observe the transition to the
bound-state contribution in Fig.6(a). However, the bound-state contribution
still shows small ocillatory behaviour shown in the small box in Fig.6(b).

Figure 6: Log plot of the numerical evaluation of the survival probability and
the pole contribution (53). Each of four plots corresponds to four points of
the parameters indicated in the phase diagram 5

####  3.2.2 Branch cut contribution

For the branch cut contribution, we resort to numerical analysis by
subtracting the analytical pole contribution shown in Eq. (53) from the full
numerical evaluation of the integral (47). Two log-log plots in the Fig.7 show
the branch cut contributions for two parameter values given in Fig. 5. In both
cases, we observe power law decay of $t^{-2}$, similar to what was observed in
the previous section.

Figure 7:

Finally, we comment on the distance between the poles and branch points. We
note that it was also analyzed in Ref.[3], in which the authors identify the
distance with a power law decay of the survival probability. We note that this
structure differs from what is seen in Section 3.1, where the distance was
between the two branch points rather than the branch point and the pole. The
distance can be associated with our phase diagram with the quantities
$e^{-|x_{1}-x_{2}|}$. This means that the lighter the green is on the phase
diagram, the wider the distance between the pole and the branch point.

##  4 Spectrual structure and Markovianity: Numerics for $L<\infty$ and $m>0$

In this section, we numerically consider the case $L<\infty$ and $m\neq 0$.
Note that the exact analysis of this case is challenging as one is required to
directly analyze the poles and branch-cut structure of the integrand in Eq.
(32), which is a transcendental equation involving a square root.

We have plotted the numerical solution of the inverse Laplace transform (32)
in Fig.8. We observe from Fig.8(b) that when the gap is open with $m=1$, the
long-time behavior remains almost unchanged, even when $L\rightarrow\infty$.
When the gap is zero, we observed a variety of dynamical behavior in different
parameter regimes shown in Fig.6. The effect of $L$ can be seen only in the
short-time regime, where in the $L\rightarrow\infty$ case, the quadratic decay
becomes linear. In conclusion, we observe that the dynamical behavior observed
in Section 3 manifest collectively in Fig.8.

Figure 8: Numerical plot of the survival probability $P(t)$ Eq. (32) with
parameters $\hbar=\omega_{0}=m=1,g=0.5$ with different cutoff values. Panel
(a) shows the log-log plot of the quantity $1-P(t)$. Panel (b) shows the log
plot of the quantity $P(t)$

##  5 Conclusion

We have considered the Markovianity and non-Markovianity of the Dirac particle
in the context of the survival probability. For our simple toy model in Eq.
(1), we have observed that the low energy structure, such as the Dirac gap
$m$, has no effect in the short time regime as shown in Fig.8(a). On the other
hand, in the long-time regime, we observed a transition from exponential decay
to power law decay as you open the Dirac gap, as shown in Table 1. We have
further performed a parameter space analysis for $m\not=0,L\rightarrow\infty$
in Fig.2 and $m=0,L<\infty$ in Fig.5. In particular, we have identified the
transition from the periodic to bound-state behavior of the survival
probability in $m=0,L<\infty$ case.

For future perspective, we aim to extend our result to the mixed state, where
Schrödinger equation is replaced with the von Neumann equation. Furthermore,
the control of long-time dynamics by the gap structure of the spectrum may
have potential applications in quantum information and quantum optics.

##  Appendix A Quantum Zeno Effect

The quantum Zeno effect is a seemingly paradoxical phenomenon in quantum
mechanics, theoretically proposed by Misra and Sudarshan [14]. It states that
a decaying quantum state under continuous observation leads to a non-decaying
probability over time. To elucidate, consider a Hamiltonian $H$ with a bounded
spectrum and an unstable state $\ket{M}$ as the initial state. Defining the
survival probability $Q(t):=|\bra{M}e^{-iHt/\hbar}\ket{M}|^{2}$. Assuming the
Hamiltonian’s expectation value with respect to the initial state $\ket{M}$ is
finite: $\bra{M}H\ket{M}<\infty$, it can be deduced that $dQ(t)/dt|_{t=0}=0$.
By segmenting the time interval into short $n$ time slices and performing a
projective measurement at each slice, the survival probability becomes
$Q(t/n)^{n}$. In the scenario of continuous measurement, as $n$ approaches
infinity, the following is derived:

$\lim_{n\rightarrow\infty}Q\left(\frac{t}{n}\right)^{n}=\lim_{n\rightarrow\infty}\left(1+\frac{1}{2}\frac{d^{2}Q(0)}{dt^{2}}\left(\frac{t}{n}\right)^{2}\right)^{n}=1.$
(56)

This suggests that continuous observation ensures the initial state $\ket{M}$
remains undecayed, given the survival probability nearing one.

##  Appendix B Contour deformation

This section shows how two contours in Fig.1 are equivalent. We begin by
adding a large contour shown as a dotted line in Fig.9(a). This large contour
does not contribute to the overall integral (33) because of the exponential in
the integrand $e^{zt/\hbar}$, which exponentially suppresses the additional
contour.

Next, we add extra contour around the branch cut, shown as solid lines in
Fig.9(b). Becasue of the extra contour, the overall integral (33) now has an
extra contribution:

$\displaystyle\int_{-iR+\sigma}^{iR+\sigma}\dots+\int_{\text{BC}}\dots~{}.$
(57)

Since the contour shown in Fig.9(b) is closed, we can use the Cauchy theorem
to deform the contour to a closed contour around the pole shown in Fig.9(d).

Finally, we must subtract the extra contribution from the branch cut by adding
another contour around the branch cut shown in Fig.9(c), but in the opposite
direction to the one added in Fig.9(b). The resulting contour shown in
Fig.9(b) is our final result.

Figure 9: Showing the process of contour deformation to show the equivalence
of two contours in Fig.1. For detail explanation of the process, see appendix
B.

##  Appendix C Laplace transform

We begin with the master equation (31). Instead of trying to find the memory
kernel directly, we perform the Laplace transform
$\text{Lap}\left[\cdot\right]$, obtaining

$\displaystyle
z\overline{\Phi}(z)-\Phi(0)=-\frac{g^{2}}{\hbar^{2}}\overline{\Phi}(z)\text{Lap}\left[\mathcal{K}(t)\right](z),$
(58)

where $\overline{\Phi}(z):=\text{Lap}\left[\Phi(t)\right](z)$.

Let us look at the Laplace transform of the memory kernel in detail:

$\displaystyle\text{Lap}\left[\mathcal{K}(t)\right](z)$ $\displaystyle=$
$\displaystyle\text{Lap}\left[e^{\frac{i}{\hbar}\omega_{0}t}\int_{-L}^{L}dp~{}\cos\left(\omega_{p}\frac{t}{\hbar}\right)\right](z)$
$\displaystyle=$
$\displaystyle\text{Lap}\left[\int_{-L}^{L}dp~{}\cos\left(\omega_{p}\frac{t}{\hbar}\right)\right]\left(z-\frac{i}{\hbar}\omega_{0}\right)$
$\displaystyle=$
$\displaystyle\int_{-L}^{L}dp~{}\text{Lap}\left[\cos\left(\omega_{p}\frac{t}{\hbar}\right)\right]\left(z-\frac{i}{\hbar}\omega_{0}\right).$

In the second equality, we have used one of the properties of the Laplace
transform, while in the last equality, we have used the definition of the
Laplace transform. The Laplace transform of the cosine function is known to be

$\displaystyle\text{Lap}\left[\cos\left(\omega_{p}\frac{t}{\hbar}\right)\right](\tilde{z})=\frac{\tilde{z}\hbar^{2}}{\tilde{z}^{2}\hbar^{2}+\omega_{p}^{2}}=\frac{\tilde{z}\hbar^{2}}{\tilde{z}^{2}\hbar^{2}+m^{2}+p^{2}},$
(59)

where $\tilde{z}\equiv z-(i/\hbar)\omega_{0}$. Integrating this quantity is
much easier than the memory kernel (31), and we thus find

$\displaystyle\int_{-L}^{L}dp~{}\frac{\tilde{z}\hbar^{2}}{\tilde{z}^{2}\hbar^{2}+m^{2}+p^{2}}=\frac{2\hbar^{2}\tilde{z}\text{Arctan}\left(\frac{L}{\sqrt{m^{2}+\tilde{z}^{2}\hbar^{2}}}\right)}{\sqrt{m^{2}+\tilde{z}^{2}\hbar^{2}}}.$
(60)

Then the Laplace transformed master equation (58) becomes

$\displaystyle
z\overline{\Phi}(z)-\Phi(0)=-\frac{g^{2}}{\hbar^{2}}\frac{2\hbar^{2}\tilde{z}\text{Arctan}\left(\frac{L}{\sqrt{m^{2}+\tilde{z}^{2}\hbar^{2}}}\right)}{\sqrt{m^{2}+\tilde{z}^{2}\hbar^{2}}}\overline{\Phi}(z).$
(61)

Rearranging the above equation, we find

$\displaystyle\overline{\Phi}(z)=\left(z+2g^{2}\frac{\tilde{z}\text{Arctan}\left(\frac{L}{\sqrt{m^{2}+\tilde{z}^{2}\hbar^{2}}}\right)}{\sqrt{m^{2}+\tilde{z}^{2}\hbar^{2}}}\right)^{-1}\Phi(0).$
(62)

The final solution is found by the inverse Laplace transform

$\displaystyle\Phi(t)=\text{Lap}^{-1}\left[\overline{\Phi}(z)\right](t).$ (63)

Once again, utilizing the shifting relation of the Laplace transform
$\text{Lap}^{-1}\left[f(z-a)\right](t)=e^{at}\text{Lap}^{-1}\left[f(z)\right](t)$,
we have

$\displaystyle\Phi(t)=e^{\frac{i}{\hbar}\omega_{0}t}\text{Lap}^{-1}\left[\overline{\Phi}\left(z+\frac{i}{\hbar}\omega_{0}\right)\right](t),$
(64)

where the function inside the inverse Laplace transform is

$\displaystyle\overline{\Phi}\left(z+\frac{i}{\hbar}\omega_{0}\right)=\left(z+i\frac{\omega_{0}}{\hbar}+2g^{2}\frac{z\text{Arctan}\left(\frac{L}{\sqrt{m^{2}+z^{2}\hbar^{2}}}\right)}{\sqrt{m^{2}+z^{2}\hbar^{2}}}\right)^{-1}\Phi(0).$
(65)

The inverse Laplace transform can be evaluated by the Bromwich integral, and
we thereby find the formal solution to the master equation (21):

$\displaystyle\Phi(t)$ $\displaystyle=$ $\displaystyle
e^{\frac{i}{\hbar}\omega_{0}t}\lim_{R\rightarrow\infty}\int_{-iR+\sigma}^{iR+\sigma}\frac{dz}{2\pi
i}e^{zt}\left(z+i\frac{\omega_{0}}{\hbar}+2g^{2}\frac{z\text{Arctan}\left(\frac{L}{\sqrt{m^{2}+z^{2}\hbar^{2}}}\right)}{\sqrt{m^{2}+z^{2}\hbar^{2}}}\right)^{-1}\Phi(0)$
(66) $\displaystyle=$ $\displaystyle
e^{\frac{i}{\hbar}\omega_{0}t}\lim_{R\rightarrow\infty}\int_{-iR+\sigma}^{iR+\sigma}\frac{dz}{2\pi
i}e^{z\frac{t}{\hbar}}\left(z+i\omega_{0}+2g^{2}\frac{z\text{Arctan}\left(\frac{L}{\sqrt{m^{2}+z^{2}}}\right)}{\sqrt{m^{2}+z^{2}}}\right)^{-1}\Phi(0)$

where in the last equality, we have resaled the integration variable by
$\hbar$.

## References

  * [1] H.-P. Breuer, F. Petruccione, et al.  The theory of open quantum systems.  Oxford University Press on Demand, 2002.
  * [2] C. Chiu, E. Sudarshan, and B. Misra.  Time evolution of unstable quantum states and a resolution of zeno’s paradox.  Physical Review D, 16(2):520, 1977.
  * [3] S. Garmon, T. Petrosky, L. Simine, and D. Segal.  Amplification of non-markovian decay due to bound state absorption into continuum.  Fortschritte der Physik, 61(2-3):261–275, 2013.
  * [4] B. Garraway.  Nonperturbative decay of an atomic system in a cavity.  Physical Review A, 55(3):2290, 1997.
  * [5] A. González-Tudela and J. I. Cirac.  Exotic quantum dynamics and purely long-range coherent interactions in dirac conelike baths.  Physical Review A, 97, 4 2018.
  * [6] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan.  Completely positive dynamical semigroups of ¡i¿n¡/i¿ -level systems.  Journal of Mathematical Physics, 17:821–825, 5 1976.
  * [7] S. Gurvitz.  Single-electron approach for time-dependent electron transport.  Physica Scripta, 2015, 10 2015.
  * [8] W. M. Itano, D. J. Heinzen, J. J. Bollinger, and D. J. Wineland.  Quantum zeno effect.  Physical Review A, 41:2295–2300, 3 1990.
  * [9] L. A. Khalfin.  Contribution to the decay theory of a quasi-stationary state, 1958.
  * [10] A. G. Kofman, G. Kurizki, and B. Sherman.  Spontaneous and induced atomic decay in photonic band structures.  Journal of Modern Optics, 41:353–384, 2 1994.
  * [11] G. Lindblad.  Mathematical physics on the generators of quantum dynamical semigroups, 1976.
  * [12] I. Mező.  On the structure of the solution set of a generalized euler–lambert equation.  Journal of Mathematical Analysis and Applications, 455(1):538–553, 2017.
  * [13] F. Minganti, A. Miranowicz, R. W. Chhajlany, and F. Nori.  Quantum exceptional points of non-hermitian hamiltonians and liouvillians: The effects of quantum jumps.  Physical Review A, 100(6):062131, 2019.
  * [14] B. Misra and E. G. Sudarshan.  The zeno’s paradox in quantum theory.  Journal of Mathematical Physics, 18(4):756–763, 1977.
  * [15] S. Nakajima.  On quantum theory of transport phenomena-steady dijfusion, 1958.
  * [16] A. Nishino and N. Hatano.  Exact time-evolving scattering states in open quantum-dot systems with an interaction: Discovery of time-evolving resonant states.  Journal of Physics A: Mathematical and Theoretical, 57(24):245302, 2024.
  * [17] R. E. A. C. Paley and N. Wiener.  Fourier transforms in the complex domain.  Monatshefte für Mathematik und Physik, 44:A8–A9, 12 1936.
  * [18] J. Redondo-Yuste, M. B. D. Paz, P. A. Huidobro, and A. González-Tudela.  Quantum electrodynamics in anisotropic and tilted dirac photonic lattices.  New Journal of Physics, 23, 10 2021.
  * [19] A. Rivas and S. F. Huelga.  Open quantum systems, volume 10.  Springer, 2012.
  * [20] C. Rothe, S. I. Hintschich, and A. P. Monkman.  Violation of the exponential-decay law at long times.  Physical Review Letters, 96:163601, 4 2006.
  * [21] E. Rutherford and F. Soddy.  Xxxiii.—the radioactivity of thorium compounds. i. an investigation of the radioactive emanation.  Journal of the Chemical Society, Transactions, 81:321–350, 1902.
  * [22] R. Zwanzig.  Ensemble method in the theory of irreversibility.  The Journal of Chemical Physics, 33:1338–1341, 1960.