# A note on the probabilistic stability
of randomized Taylor schemes

Tomasz Bochacik  AGH University of Science and Technology, Faculty of Applied
Mathematics,
Al. A. Mickiewicza 30, 30-059 Kraków, Poland<EMAIL_ADDRESS>
###### Abstract.

We study the stability of randomized Taylor schemes for ODEs. We consider
three notions of probabilistic stability: asymptotic stability, mean-square
stability, and stability in probability. We prove fundamental properties of
the probabilistic stability regions and benchmark them against the absolute
stability regions for deterministic Taylor schemes.

Key words: randomized Taylor schemes, mean-square stability, asymptotic
stability, stability in probability

MSC 2010: 65C05, 65L05, 65L20

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
  3. 3 Main results
  4. 4 Conclusions
  5. A Auxiliary lemmas

##  1\. Introduction

The study of randomized algorithms approximating the solutions of initial
value problems for ODEs dates back to early 1990s, cf. [16, 17].

So far, the main focus has been on convergence of randomized algorithms, see
for example [6, 10, 14]. Randomized algorithms tend to converge faster than
their deterministic counterparts, especially for the problems of low
regularity. Error bounds are usually established using certain martingale
inequalities and classical tools such as Gronwall’s inequality. In many
papers, error analysis was combined with the discussion of algorithms’
optimality (in the Information Based Complexity sense), cf. [3, 4, 7, 8, 11,
12, 13], also in the setting of inexact information, cf. [1, 2]. Randomized
Taylor schemes, which will be of particular interest in this paper, were shown
in [8] to achieve the optimal rate of convergence under mild regularity
conditions.

Other aspects of randomized algorithms, such as stability, have been largely
omitted. The aim of this paper is to make a step towards filling this gap.

The stability of deterministic algorithms for ODEs has been comprehensively
studied in the literature, see, for example, [5]. In the stability analysis,
we consider a test problem which is simple enough but retains features present
in a wider class of problems. Then, we investigate for which choices of the
step-size, the method reproduces the characteristics of the test equation, cf.
[9]. In the context of ODEs, we usually take a linear, scalar, and autonomous
test problem.

The same test problem is used for stability analysis of randomized algorithms
for ODEs. However, the approximated solution generated by a randomized method
is random. Hence, analysis of its behaviour in infinity depends on the type of
convergence. This naturally leads to notions of the mean-square stability and
the asymptotic (almost-sure) stability, which have been previously considered
in [9, 15] in the context of stochastic differential equations. This
framework, enriched with the notion of stability in probability, was used in
[1, 2] to characterize stability regions of randomized Euler schemes and the
randomized two-stage Runge-Kutta scheme.

This paper, according to our best knowledge, is the first attempt to apply the
concept of probabilistic stability to higher-order randomized methods for
ODEs, namely to the family of randomized Taylor schemes defined in [8]. Since
these methods do not involve implicitness, they will not be A-stable. However,
we can characterize probabilistic stability regions of these methods in a
quite detailed way. We establish their basic properties such as openness,
boundedness, symmetry. Moreover, we study inclusions between them and compare
them to the reference sets corresponding to deterministic methods. Finally, we
provide counterexamples for some hypothetical properties which do not hold for
probabilistic stability regions of randomized Taylor schemes in general (that
is, for the method of any order).

In Section 2 we give basic definitions and introduce the notation. In
particular, we recall the definitions of the family of randomized Taylor
schemes and of the probabilistic stability regions. In Section 3, which is the
main part of this paper, we characterize probabilistic stability regions for
randomized Taylor schemes. Conclusions are discussed in Section 4. Finally, in
Appendix A we prove some technical lemmas.

##  2\. Preliminaries

###  2.1. The family of randomized Taylor schemes

We deal with initial value problems of the following form:

$\left\\{\begin{array}[]{ll}z^{\prime}(t)=f(t,z(t)),\ t\in[a,b],\\\
z(a)=\eta,\end{array}\right.$ (1)

where $-\infty<a<b<\infty$, $\eta\in\mathbb{R}^{d}$,
$f\colon[a,b]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$, $d\in\mathbb{Z}_{+}$.

We use the definition of the family of randomized Taylor schemes given in [8].
We fix $n\in\mathbb{N}$, $n\geq 2$ and put $h=\frac{b-a}{n}$, $t_{j}=a+jh$ for
$j\in\\{0,1,\ldots,n\\}$, $\theta_{j}=t_{j-1}+\tau_{j}h$, $\tau_{j}\sim
U(0,1)$ for $j\in\\{1,\ldots,n\\}$. We assume that the family of random
variables $\\{\tau_{1},\ldots,\tau_{n}\\}$ is independent.

Let $r\in\mathbb{N}$ (we assume that $0\in\mathbb{N}$). We set
$v^{r}_{0}=\eta$. If $k\in\\{1,\ldots,n\\}$ and $v^{r}_{k-1}$ is already
defined, we consider the following local problem:

$\left\\{\begin{array}[]{ll}(u^{r}_{k})^{\prime}(t)=f(t,u^{r}_{k}(t)),\
t\in[t_{k-1},t_{k}],\\\ u^{r}_{k}(t_{k-1})=v^{r}_{k-1}.\end{array}\right.$ (2)

We define

$p^{r}_{k}(t)=\sum_{j=0}^{r+1}\frac{(u^{r}_{k})^{(j)}(t_{k-1})}{j!}(t-t_{k-1})^{j}\mathbbm{1}_{[t_{k-1},t_{k}]}(t)$
(3)

and

$v^{r}_{k}=p^{r}_{k}(t_{k})+h\cdot\bigl{(}f(\theta_{k},p^{r}_{k}(\theta_{k}))-(p^{r}_{k})^{\prime}(\theta_{k})\bigr{)}.$
(4)

Note that for $t\in(t_{k-1},t_{k})$, where $k\in\\{1,\ldots,n\\}$, we have

$(p^{r}_{k})^{\prime}(t)=\sum_{j=0}^{r+1}\frac{(u^{r}_{k})^{(j)}(t_{k-1})}{(j-1)!}(t-t_{k-1})^{j-1}.$
(5)

The algorithm returns the sequence $\bigl{(}v_{k}^{r}\bigr{)}_{k=0}^{n}$,
which approximates values of the exact solution $z$ of (1) at points
$t_{0},\ldots,t_{n}$. If we neglect the second term in the right-hand side of
(4), we get the classical deterministic Taylor scheme.

Note that for $r=0$ we obtain the definition of randomized two-stage Runge-
Kutta scheme given in [1, 14]. Stability analysis for this algorithm was
performed in [1].

###  2.2. Probabilistic stability of randomized Taylor schemes

Let us consider the classical test problem

$\left\\{\begin{array}[]{ll}z^{\prime}(t)=\lambda z(t),\ t\geq 0,\\\
z(0)=\eta\end{array}\right.$ (6)

with $\lambda\in\mathbb{C}$ and $\eta\neq 0$. The exact solution of (6) is
$z(t)=\eta\exp(\lambda t)$. We note that

$\lim\limits_{t\to\infty}z(t)=0\ \hbox{iff}\ \Re(\lambda)<0.$

For a fixed step-size $h>0$, we apply the scheme (4) with the mesh $t_{k}=kh$,
$k\in\mathbb{N}$, to the test problem (6). The local problem (2) takes on the
following form:

$\left\\{\begin{array}[]{ll}(u^{r}_{k})^{\prime}(t)=\lambda u^{r}_{k}(t),\
t\in[t_{k-1},t_{k}],\\\ u^{r}_{k}(t_{k-1})=v^{r}_{k-1}.\end{array}\right.$ (7)

As a result, we obtain a sequence $(v^{r}_{k})_{k=0}^{\infty}$, which
approximates the values of $z$ at $t_{k}$ ($k\in\mathbb{N}$) and whose values
are given by

$v^{r}_{k}=\eta\cdot\prod_{l=1}^{k}f_{r,\lambda h}(\tau_{l}),$ (8)

where $\tau_{1},\tau_{2},\ldots$ are independent random variables with uniform
distribution on $[0,1]$, and

$f_{r,z}\colon\mathbb{R}\ni
t\mapsto\sum_{j=0}^{r+1}\frac{z^{j}}{j!}+\frac{t^{r+1}}{(r+1)!}z^{r+2}\in\mathbb{C}.$
(9)

In fact, using (7) and proceeding by induction with respect to $j$, we get

$(u^{r}_{k})^{(j)}(t_{k-1})=\lambda^{j}v^{r}_{k-1}.$

Thus, by taking $t=t_{k}$ and $t=\theta_{k}$ in (3) and (5), we get

$\displaystyle p^{r}_{k}(t_{k})$ $\displaystyle=\sum_{j=0}^{r+1}\frac{(\lambda
h)^{j}}{j!}v^{r}_{k-1},$ $\displaystyle p^{r}_{k}(\theta_{k})$
$\displaystyle=\sum_{j=0}^{r+1}\frac{(\lambda h\tau_{k})^{j}}{j!}v^{r}_{k-1},$
$\displaystyle(p^{r}_{k})^{\prime}(\theta_{k})$
$\displaystyle=\sum_{j=1}^{r+1}\frac{\lambda^{j}(h\tau_{k})^{j-1}}{(j-1)!}v^{r}_{k-1}=\lambda\sum_{j=0}^{r}\frac{(\lambda
h\tau_{k})^{j}}{j!}v^{r}_{k-1}.$

By (4) and the above three lines, we obtain the following recurrence:

$\displaystyle v^{r}_{k}$ $\displaystyle=p^{r}_{k}(t_{k})+\lambda
hp^{r}_{k}(\theta_{k})-h(p^{r}_{k})^{\prime}(\theta_{k})=\Bigl{(}\sum_{j=0}^{r+1}\frac{(\lambda
h)^{j}}{j!}+\lambda h\cdot\frac{(\lambda
h\tau_{k})^{r+1}}{(r+1)!}\Bigr{)}\cdot v^{r}_{k-1},$

which leads to (8).

Similarly as in [1], we consider three sets

$\displaystyle\mathcal{R}^{r}_{MS}=\\{\lambda h\in\mathbb{C}\colon
v^{r}_{k}\to 0\ \hbox{in $L^{2}(\Omega)$ as}\ k\to\infty\\},$ (10)
$\displaystyle\mathcal{R}^{r}_{AS}=\\{\lambda h\in\mathbb{C}\colon
v^{r}_{k}\to 0\ \hbox{almost surely as}\ k\to\infty\\},$ (11)
$\displaystyle\mathcal{R}^{r}_{SP}=\left\\{\lambda h\in\mathbb{C}\colon
v^{r}_{k}\to 0\text{ in probability as}\ k\to\infty\right\\},$ (12)

where we call $\mathcal{R}^{r}_{MS}$ the region of mean-square stability,
$\mathcal{R}^{r}_{AS}$ – the region of asymptotic stability, and
$\mathcal{R}^{r}_{SP}$ – the region of stability in probability.

######  Remark 1.

If $h\lambda\in\mathcal{R}^{r}_{AS}$, the approximated solution of the test
problem converges to $0$ for virtually every fixing of
$\tau_{1},\tau_{2},\ldots$ from the interval $[0,1]$. That is, for almost
every $\omega\in\Omega$ and for every $\varepsilon>0$, there exists
$k_{0}=k_{0}(\omega,\varepsilon)$ such that $|v_{k}^{r}(\omega)|<\varepsilon$
for every $k\geq k_{0}$.

The mean-square stability implies that for each pre-specfied threshold
$\varepsilon>0$ and each pre-specified probability level $\delta\in(0,1)$, the
point-wise approximations $v^{r}_{k}$ (for sufficiently big $k$) are bounded
by $\varepsilon$ with probability at least $1-\delta$. In fact, let us assume
that $h\lambda\in\mathcal{R}^{r}_{MS}$ and let us set $\varepsilon>0$ and
$\delta\in(0,1)$. Then there exists
$k^{\prime}_{0}=k^{\prime}_{0}(\varepsilon,\delta)\in\mathbb{Z}_{+}$ such that
for each $k\geq k^{\prime}_{0}$,
$\|v^{r}_{k}\|_{L^{2}(\Omega)}<\sqrt{\delta}\varepsilon$. By Markov’s
inequality,

$\mathbb{P}\bigl{(}|v^{r}_{k}|\geq\varepsilon\bigr{)}\leq\frac{\|v^{r}_{k}\|^{2}_{L^{2}(\Omega)}}{\varepsilon^{2}}<\delta.$

Hence, the asymptotic stability guarantees that for each specific run of the
algorithm, the approximated solution will converge to the exact solution in
infinity (provided that $h\lambda\in\mathcal{R}^{r}_{AS}$). However, the pace
of convergence may significantly vary for different fixings of
$\tau_{1},\tau_{2},\ldots$. On the other hand, the mean-square stability
provides an insight on whether we can expect consistent asymptotic behaviour
in many independent runs of the algorithm.

##  3\. Main results

For fixed $r\in\mathbb{N}$, let

$\displaystyle F_{r}\colon$ $\displaystyle\mathbb{C}\ni
z\mapsto\mathbb{E}|f_{r,z}(\tau)|^{2}=\int\limits_{0}^{1}|f_{r,z}(t)|^{2}\,\mathrm{d}t$
(13) $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \
=\int\limits_{0}^{1}\Bigl{|}\sum_{j=0}^{r+1}\frac{z^{j}}{j!}+\frac{t^{r+1}}{(r+1)!}z^{r+2}\Bigr{|}\,\mathrm{d}t\in[0,\infty),$
$\displaystyle G_{r}\colon$ $\displaystyle\mathbb{C}\ni
z\mapsto\mathbb{E}\bigl{(}\ln|f_{r,z}(\tau)|\bigr{)}=\int\limits_{0}^{1}\ln|f_{r,z}(t)|\,\mathrm{d}t\in\mathbb{R},$
(14)

where $\tau\sim U([0,1])$ and $f_{r,z}$ is given by (9). From (8), (10) and
(13) we obtain

$\mathcal{R}_{MS}^{r}=\bigl{\\{}z\in\mathbb{C}\colon
F_{r}(z)<1\bigr{\\}}=\left\\{z\in\mathbb{C}\colon\int\limits_{0}^{1}\Bigl{|}\sum_{j=0}^{r+1}\frac{z^{j}}{j!}+\frac{t^{r+1}}{(r+1)!}z^{r+2}\Bigr{|}\,\mathrm{d}t<1\right\\}.$
(15)

For further analysis of the probabilistic regions of randomized Taylor
schemes, we need Lemmas 1–4. They are formulated and proven in Appendix A.

By Lemma 1 and Corollary 1 in [1], regions $\mathcal{R}^{r}_{AS}$ and
$\mathcal{R}^{r}_{SP}$ are equal and can be expressed as

$\mathcal{R}^{r}_{AS}=\mathcal{R}^{r}_{SP}=\bigl{\\{}z\in\mathbb{C}\colon
G_{r}(z)<0\bigr{\\}}.$ (16)

Apart from $\mathcal{R}_{MS}^{r},\mathcal{R}_{AS}^{r}$ and
$\mathcal{R}_{SP}^{r}$, we consider the following reference set:

$\mathcal{R}_{ref}^{r}=\Bigl{\\{}z\in\mathbb{C}\colon\eta\prod_{l=1}^{k}\mathbb{E}f_{r,z}(\tau_{l})\to
0\text{ as
}k\to\infty\Bigr{\\}}=\Bigl{\\{}z\in\mathbb{C}\colon\Bigl{|}\sum_{j=0}^{r+2}\frac{z^{j}}{j!}\Bigr{|}<1\Bigr{\\}}.$
(17)

Note that the reference set $\mathcal{R}_{ref}^{r}$ is the absolute stability
region for the deterministic Taylor scheme of order $r+2$. Hence, our analysis
is consistent with [1], where probabilistic stability regions of the
randomized two-stage Runge-Kutta scheme (i.e., the randomized Taylor scheme
with $r=0$) were benchmarked against the absolute stability region of the mid-
point method (i.e., the deterministic Taylor scheme with $r=1$).

Now we are ready to establish the main results of this paper – Theorem 1 and
Theorem 2. They extend Theorem 3 and Theorem 4 from [1], which covered only
the case of $r=0$, onto the case of any $r\in\mathbb{N}$. However, as of now
we have not managed to generalize the results from [1] related to stability
intervals.

######  Theorem 1.

For each $r\in\mathbb{N}$, the sets $\mathcal{R}^{r}_{MS}$,
$\mathcal{R}^{r}_{AS}$ and $\mathcal{R}^{r}_{ref}$ are open and symmetric with
respect to the real axis.

###### Proof.

Sets $\mathcal{R}^{r}_{MS}=F_{r}^{-1}((-\infty,1))$ and
$\mathcal{R}^{r}_{AS}=G_{r}^{-1}((-\infty,0))$ are open due to Lemma 2 and
Lemma 4, respectively. Since the function

$\hat{F}_{r}\colon\mathbb{C}\ni
z\mapsto\Bigl{|}\sum_{j=0}^{r+2}\frac{z^{j}}{j!}\Bigr{|}\in[0,\infty)$

is continuous, set $\mathcal{R}^{r}_{ref}=\hat{F}_{r}^{-1}((-\infty,1))$ is
open as well.

Note that $|f_{r,z}(t)|=|f_{r,\bar{z}}(t)|$ for all $z\in\mathbb{C}$ and
$t\in[0,1]$, which implies that $F_{r}(z)=F_{r}(\bar{z})$ and
$G_{r}(z)=G_{r}(\bar{z})$ for all $z\in\mathbb{C}$, cf. (9), (13) and (14).
This combined with (15) and (16) immediately gives the symmetry of
$\mathcal{R}^{r}_{MS}$ and $\mathcal{R}^{r}_{AS}$ with respect to the real
axis. The same property for $\mathcal{R}^{r}_{ref}$ follows from the fact that

$\Bigl{|}\sum_{j=0}^{r+2}\frac{z^{j}}{j!}\Bigr{|}=\Bigl{|}\sum_{j=0}^{r+2}\frac{\bar{z}^{j}}{j!}\Bigr{|}$

for all $z\in\mathbb{C}$. ∎

The following Theorem 2 shows that the mean-square stability is a stronger
property than the asymptotic stability. Furthermore, probabilistic stability
regions of the randomized Taylor scheme for any $r\in\mathbb{N}$ are bounded.
This implies that none of randomized Taylor schemes is A-stable in any of the
considered probabilistic senses (i.e., the left complex half-plane is not
contained in any of the considered stability regions). However, the left half-
plane is contained in the sum over $r$ of stability regions.

######  Theorem 2.

For each $r\in\mathbb{N}$, there exists $\gamma_{r}\in(0,\infty)$ such that

$\mathcal{R}^{r}_{MS}\subset\mathcal{R}^{r}_{ref}\cap\mathcal{R}^{r}_{AS}\subset\mathcal{R}^{r}_{ref}\cup\mathcal{R}^{r}_{AS}\subset\\{z\in\mathbb{C}\colon|z|<\gamma_{r}\\}.$
(18)

Moreover,

$\mathbb{C}_{-}\subset\bigcup_{r=0}^{\infty}\mathcal{R}^{r}_{MS}.$ (19)

###### Proof.

Inclusion $\mathcal{R}_{MS}^{r}\subset\mathcal{R}_{ref}^{r}$ for all
$r\in\mathbb{N}$ follows from the following inequality:

$F_{r}(z)=\mathbb{E}\bigl{|}f_{r,z}(\tau)\bigr{|}^{2}\geq\bigl{|}\mathbb{E}f_{r,z}(\tau)\bigr{|}^{2}=\Bigl{|}\sum_{j=0}^{r+1}\frac{z^{j}}{j!}+\frac{\mathbb{E}\tau^{r+1}}{(r+1)!}z^{r+2}\Bigr{|}^{2}=\Bigl{|}\sum_{j=0}^{r+2}\frac{z^{j}}{j!}\Bigr{|}^{2}.$

We used the fact that $\mathbb{E}|Z|^{2}\geq|\mathbb{E}Z|^{2}$ for any complex
random variable $Z$. Since convergence in $L^{2}(\Omega)$ implies convergence
in probability, we have
$\mathcal{R}_{MS}^{r}\subset\mathcal{R}_{SP}^{r}=\mathcal{R}_{AS}^{r}$ for all
$r\in\mathbb{N}$, cf. (10), (12) and (16). As a result, we obtain the first
inclusion in (18).

Region $\mathcal{R}_{ref}^{r}$ is bounded for each $r\in\mathbb{N}$ because

$\Bigl{|}\sum_{j=0}^{r+2}\frac{z^{j}}{j!}\Bigr{|}\geq\frac{|z|^{r+2}}{(r+2)!}-\sum_{j=0}^{r+1}\frac{|z|^{j}}{j!}$

and the right-hand side of the above inequality tends to infinity when
$|z|\to\infty$. Thus, there exists $\gamma^{1}_{r}>0$ such that
$\Bigl{|}\sum_{j=0}^{r+2}\frac{z^{j}}{j!}\Bigr{|}\geq 1$ for all
$z\in\mathbb{C}$ with $|z|\geq\gamma^{1}_{r}$.

To show that $\mathcal{R}_{AS}^{r}$ is bounded, let us express $G_{r}(z)$ for
$z\in\mathbb{C}\setminus\\{0\\}$ in the same fashion as in the proof of Lemma
4:

$G_{r}(z)=\ln\Bigl{(}\frac{|z|^{r+2}}{(r+1)!}\Bigr{)}+H_{r}(h_{r}(z))\geq(r+2)\ln|z|-\ln((r+1)!)+\inf_{z\in\mathbb{C}}H_{r}(z),$
(20)

cf. (28), (29) and (30). Note that
$\displaystyle\inf_{z\in\mathbb{C}}H_{r}(z)$ is finite because

$H_{r}(z)\geq\int\limits_{0}^{1}\ln\bigl{|}|z|-t^{r+1}\bigr{|}\,\mathrm{d}t\geq\int\limits_{0}^{1}\ln(2-t^{r+1})\,\mathrm{d}t\geq
0$

for $|z|>2$ and $\displaystyle\inf_{|z|\leq 2}H_{r}(z)$ is finite by Lemma 4
and the Weierstrass extreme value theorem. The right-hand side of (20) tends
to infinity when $|z|\to\infty$. Hence, there exists $\gamma^{2}_{r}>0$ such
that $G_{r}(z)\geq 0$ for all $z\in\mathbb{C}$ with $|z|\geq\gamma^{2}_{r}$.
Taking $\gamma_{r}=\max\\{\gamma^{1}_{r},\gamma^{2}_{r}\\}$ leads to the third
inclusion in (18).

To see (19), let us consider $z\in\mathbb{C}$ such that
$z\notin\bigcup_{r=0}^{\infty}\mathcal{R}^{r}_{MS}$. Then $F_{r}(z)\geq 1$ for
all $r\in\mathbb{N}$ and as a result

$1\leq\limsup_{r\to\infty}\int\limits_{0}^{1}\bigl{|}f_{r,z}(t)\bigr{|}^{2}\,\mathrm{d}t\leq\int\limits_{0}^{1}\limsup_{r\to\infty}\bigl{|}f_{r,z}(t)\bigr{|}^{2}\,\mathrm{d}t=\int\limits_{0}^{1}|e^{z}|^{2}\,\mathrm{d}t=e^{2\Re(z)}.$
(21)

Hence, $\Re(z)\geq 0$ and (19) follows. The second inequality in (21) is based
on Fatou’s lemma. Note that $f_{r,z}$ is continuous in $\mathbb{R}$, which
guarantees that the function $[0,1]\ni
t\mapsto\bigl{|}f_{r,z}(t)\bigr{|}^{2}\in[0,\infty)$ is continuous as well and
thus Borel measurable. In the third passage in (21) we use the fact that
$\displaystyle\lim_{r\to\infty}f_{r,z}(t)$ exists and is equal to $e^{z}$ for
each $t\in[0,1]$. ∎

In Figure 1, we plot $\mathcal{R}^{r}_{ref}$, $\mathcal{R}^{r}_{MS}$, and
$\mathcal{R}^{r}_{AS}$ for $r\in\\{0,1,2,3,4\\}$. Based on these plots and
some calculations, we have rejected several hypotheses about potential
properties of these regions. Counterexamples are provided in Remarks 2–6.

(a) Contours of $\mathcal{R}^{r}_{ref}$.

(b) Contours of $\mathcal{R}^{r}_{MS}$.

(c) Contours of $\mathcal{R}^{r}_{AS}$.

Figure 1. Contours of reference sets, regions of mean-sqare stability, and
regions of asymptotic stability of the randomized Taylor schemes for $r=0$
(red line), $r=1$ (purple line), $r=2$ (orange line), $r=3$ (blue line), and
$r=4$ (gray line).

######  Remark 2.

None of the following inclusions holds in general (for every
$r\in\mathbb{N}$):

  * a)

$\mathcal{R}^{r}_{ref}\subset\mathcal{R}^{r+1}_{ref}$;

  * b)

$\mathcal{R}^{r}_{MS}\subset\mathcal{R}^{r+1}_{MS}$;

  * c)

$\mathcal{R}^{r}_{AS}\subset\mathcal{R}^{r+1}_{AS}$.

Let us consider $z_{a}=-0.6+2.8\,i$, $z_{b}=-0.03+1.9\,i$, and
$z_{c}=-0.25+2.75\,i$. First two of these points are represented as the
intersection of dashed lines in Figure 2. We have

$\Bigl{|}\sum_{j=0}^{4}\frac{z_{a}^{j}}{j!}\Bigr{|}^{2}=\frac{253\,409}{360\,000}<1\
\text{ and }\
\Bigl{|}\sum_{j=0}^{5}\frac{z_{a}^{j}}{j!}\Bigr{|}^{2}=\frac{5\,828\,357}{5\,625\,000}>1.$

Thus, by (17), $z_{a}\in\mathcal{R}_{ref}^{2}\setminus\mathcal{R}_{ref}^{3}$.
Furthermore,

$\displaystyle F_{2}(z_{b})$
$\displaystyle=\frac{2\,460\,549\,996\,776\,228\,711}{2\,520\,000\,000\,000\,000\,000}<1,$
$\displaystyle F_{3}(z_{b})$
$\displaystyle=\frac{531\,703\,423\,127\,449\,318\,399\,669}{518\,400\,000\,000\,000\,000\,000\,000}>1,$

which means that $z_{b}\in\mathcal{R}_{MS}^{2}\setminus\mathcal{R}_{MS}^{3}$,
cf. (15). Finally, using the scipy.integrate.quad function in Python, we
obtain the following estimates:

$G_{3}(z_{c})\approx-0.41731<0\ \ \text{ and }\ \ G_{4}(z_{c})\approx
0.06505>0.$

Thus, by (16), $z_{c}\in\mathcal{R}^{3}_{AS}\setminus\mathcal{R}^{4}_{AS}$.

(a) Contours of $\mathcal{R}^{2}_{ref}$ and $\mathcal{R}^{3}_{ref}$.

(b) Contours of $\mathcal{R}^{2}_{MS}$ and $\mathcal{R}^{3}_{MS}$.

Figure 2. Fragments of contours of reference sets and regions of mean-sqare
stability of the randomized Taylor schemes for $r=2$ (orange line) and $r=3$
(blue line).

######  Remark 3.

In general, there is no inclusion between $\mathcal{R}^{r}_{AS}$ and
$\mathcal{R}^{r}_{ref}$.

To see this, let us consider $r=0$ and take $z_{1}=-2.1$, $z_{2}=-1+1.6\,i$.
Then

$\displaystyle\Bigl{|}1+z_{1}+\frac{z_{1}^{2}}{2}\Bigr{|}=1.105>1\ \text{ and
}\ \Bigl{|}1+z_{2}+\frac{z_{2}^{2}}{2}\Bigr{|}=0.78<1.$

On the other hand, we have estimated $G_{0}(z_{1})\approx-0.07784$ and
$G_{0}(z_{2})\approx 0.13565$. By (14), (16) and (17), we obtain
$z_{1}\in\mathcal{R}^{0}_{AS}\setminus\mathcal{R}^{0}_{ref}$ and
$z_{2}\in\mathcal{R}^{0}_{ref}\setminus\mathcal{R}^{0}_{AS}$.

######  Remark 4.

In general, regions $\mathcal{R}^{r}_{ref}$, $\mathcal{R}^{r}_{MS}$ and
$\mathcal{R}^{r}_{AS}$ are not included in $\mathbb{C}_{-}$.

This inclusion is true for $r=0$, cf. Theorem 3(ii) and Theorem 4(iii) in [1],
but does not hold for $r=1$, as shown in the following example:

$F_{1}(0.01+i)=\frac{19\,772\,000\,147\,001}{20\,000\,000\,000\,000}<1,$

which combined with (15) and (18) implies that
$0.01+i\in\mathcal{R}^{1}_{MS}\cap\mathcal{R}^{1}_{ref}\cap\mathcal{R}^{1}_{AS}\cap\mathbb{C}_{+}$.

######  Remark 5.

In general, sets $\mathcal{R}^{r}_{ref}$, $\mathcal{R}^{r}_{MS}$ and
$\mathcal{R}^{r}_{AS}$ are not convex.

Let $z_{1}=0.01+i$ and $z_{2}=0.01-i$. We know that
$z_{1},z_{2}\in\mathcal{R}^{1}_{MS}\cap\mathcal{R}^{1}_{ref}\cap\mathcal{R}^{1}_{AS}$,
cf. Remark 4 and Theorem 1. On the other hand,
$(\mathcal{R}^{1}_{MS}\cup\mathcal{R}^{1}_{ref}\cup\mathcal{R}^{1}_{AS})\cap\mathbb{R}\subset(-\infty,0)$
because $f_{r,z}(t)\geq 1$ for all $r\in\mathbb{N}$, $z\in[0,\infty)$,
$t\in[0,1]$, and we have (14), (16), (17), (18). Thus,
$\frac{z_{1}+z_{2}}{2}=0.01\notin\mathcal{R}^{1}_{MS}\cup\mathcal{R}^{1}_{ref}\cup\mathcal{R}^{1}_{AS}$.

######  Remark 6.

In general, sets $\mathcal{R}^{r}_{ref}$ and $\mathcal{R}^{r}_{AS}$ are not
connected.

Let us take $z_{1}=0.75+3.5\,i$ and $z_{2}=-0.25+2.5\,i$. We have

$\displaystyle\Bigl{|}\sum_{j=0}^{6}\frac{z_{1}^{j}}{j!}\Bigr{|}^{2}$
$\displaystyle=\frac{27\,473\,196\,877\,335\,817\,540\,321}{121\,029\,087\,867\,608\,368\,152\,576}<1,$
$\displaystyle\Bigl{|}\sum_{j=0}^{6}\frac{z_{2}^{j}}{j!}\Bigr{|}^{2}$
$\displaystyle=\frac{48\,715\,333\,577\,673\,689\,545\,536\,241}{75\,643\,179\,917\,255\,230\,095\,360\,000}<1$

but

$\Bigl{|}\sum_{j=0}^{6}\frac{\Bigl{(}\frac{z_{1}+z_{2}}{2}\Bigr{)}^{j}}{j!}\Bigr{|}^{2}=\frac{9\,427\,129\,581\,150\,440\,422\,815\,049}{3\,025\,727\,196\,690\,209\,203\,814\,400}>1.$

Thus, $z_{3},z_{4}\in\mathcal{R}^{4}_{ref}$ but
$\frac{z_{3}+z_{4}}{2}\notin\mathcal{R}^{4}_{ref}$, cf. (17), which means that
the set $\mathcal{R}^{4}_{ref}$ is disconnected.

For $z_{3}=-0.5+2\,i$ and $z_{4}=0.25+3.25\,i$ we get

$G_{4}(z_{3})\approx-0.50028<0,\ \ G_{4}(z_{4})\approx-0.47024<0,\ \text{ and
}\ G_{4}\Bigl{(}\frac{z_{3}+z_{4}}{2}\Bigr{)}\approx 0.03656>0,$

which implies that $z_{3},z_{4}\in\mathcal{R}^{4}_{AS}$ but
$\frac{z_{3}+z_{4}}{2}\notin\mathcal{R}^{4}_{AS}$, see (16). Hence, the set
$\mathcal{R}^{4}_{AS}$ is disconnected.

In general, disconnectivity of the stability region would indicate that the
method’s behaviour for stiff problems is in a sense unpredictable – taking
smaller $h$ would not necessarily improve the method’s performance. However,
disconnected parts of the reference and asymptotic stability regions are
observed in the right complex half-plane, which is out of interest in the
context of A-stability. We conjecture that the intersection of each of the
aforementioned stability regions with the left half-plane is connected.

##  4\. Conclusions

We have established fundamental properties of probabilistic stability regions
for randomized Taylor schemes. In particular, we have shown that notions of
asymptotic stability and stability in probability are equivalent for this
family of schemes, cf. (16). Furthermore, we have proven openness and symmetry
of all considered stability regions (see Theorem 1), as well as their
boundedness, cf. (18) in Theorem 2.

Although randomized Taylor schemes are not A-stable for any $r\in\mathbb{N}$
and in any probabilistic sense, the union of (asymptotic or mean-square)
stability regions over all $r\in\mathbb{N}$ covers the entire left complex
half-plane, see (19) in Theorem 2. Hence, if the right-hand side function $f$
is sufficiently regular (in the most optimistic scenario, analytical), one may
increase $r$ in order to prevent rapid variation in the approximated solution
of a stiff problem. Otherwise, applicability of the methods studied in this
paper is in practice limited to non-stiff problems.

Finally, we have ruled out a number of hypotheses concerning the potential
properties of stability regions for randomized Taylor schemes, including their
monotonicity (with respect to $r$) and some other potential inclusions,
convexity, and connectivity, cf. Remarks 2–6.

## Acknowledgments

I would like to thank Professor Paweł Przybyłowicz for many inspiring
discussions while preparing this manuscript.

I am also very grateful to the anonymous reviewer for valuable comments and
suggestions that allowed me to improve the quality of the paper.

This research was funded in whole or in part by the National Science Centre,
Poland, under project 2021/41/N/ST1/00135. For the purpose of Open Access, the
author has applied a CC-BY public copyright licence to any Author Accepted
Manuscript (AAM) version arising from this submission.

##  Appendix A Auxiliary lemmas

In this section, we prove technical lemmas which are necessary to establish
equality (16), Theorem 1, and Theorem 2.

######  Lemma 1.

For each $z\in\mathbb{C}$ and each $r\in\mathbb{N}$, the random variable
$\ln|f_{r,z}(\tau)|$ is square-integrable.

###### Proof.

We will show more, i.e.,

$\int\limits_{0}^{1}\ln^{2}\bigl{|}z_{1}t^{k}+z_{2}\bigl{|}^{2}\,\mathrm{d}t<\infty$
(22)

for all $k\in\mathbb{Z}_{+}$ and $z_{1},z_{2}\in\mathbb{C}$ such that
$(z_{1},z_{2})\neq(0,0)$. Since the case $z_{1}=0$ is trivial, in the
following we consider any $z_{1},z_{2}\in\mathbb{C}$ such that $z_{1}\neq 0$.
If $-\frac{z_{2}}{z_{1}}\notin[0,1]$, then $z_{1}t^{k}+z_{2}\neq 0$ for all
$t\in[0,1]$ and (22) follows because the integrand is a continuous function of
$t$ on the interval $[0,1]$. From this point we assume that
$\alpha=-\frac{z_{2}}{z_{1}}\in[0,1]$. Then

$\int\limits_{0}^{1}\ln^{2}\bigl{|}z_{1}t^{k}+z_{2}\bigl{|}\,\mathrm{d}t=\int\limits_{0}^{1}\Bigl{(}\ln|z_{1}|+\ln\bigl{|}t^{k}-\alpha\bigl{|}\Bigr{)}^{2}\,\mathrm{d}t.$

If $\alpha=0$, both integrals

$\int\limits_{0}^{1}\ln\bigl{|}t^{k}-\alpha\bigl{|}\,\mathrm{d}t=k\int\limits_{0}^{1}\ln
t\,\mathrm{d}t\ \text{ and }\
\int\limits_{0}^{1}\bigl{(}\ln\bigl{|}t^{k}-\alpha\bigl{|}\bigr{)}^{2}\,\mathrm{d}t=k^{2}\int\limits_{0}^{1}(\ln
t)^{2}\,\mathrm{d}t$

are finite. In case of $\alpha\in(0,1]$, there exists $\beta\in(0,1]$ such
that $\alpha=\beta^{k}$ and we can write that

$\int\limits_{0}^{1}\ln\bigl{|}t^{k}-\alpha\bigl{|}\,\mathrm{d}t=\int\limits_{0}^{1}\ln|t-\beta|\,\mathrm{d}t+\int\limits_{0}^{1}\ln\Bigl{(}\beta^{k-1}+\sum_{j=0}^{k-2}t^{k-1-j}\beta^{j}\Bigr{)}\,\mathrm{d}t.$

The first integral above has one singularity but this is a well-known fact
that it is finite. For the second one, we note that the integrand is a
continuous function for $t\in[0,1]$. Furthermore, since $(a+b)^{2}\leq
2a^{2}+2b^{2}$ for any real numbers $a,b$, we obtain

$\int\limits_{0}^{1}\bigl{(}\ln\bigl{|}t^{k}-\alpha\bigl{|}\bigr{)}^{2}\,\mathrm{d}t\leq
2\int\limits_{0}^{1}\bigl{(}\ln|t-\beta|\bigr{)}^{2}\,\mathrm{d}t+2\int\limits_{0}^{1}\Bigl{[}\ln\Bigl{(}\beta^{k-1}+\sum_{j=0}^{k-2}t^{k-1-j}\beta^{j}\Bigr{)}\Bigr{]}^{2}\,\mathrm{d}t$

and we may use similar arguments as before to justify that the above integrals
are finite. ∎

######  Lemma 2.

For each $r\in\mathbb{N}$, the function $F_{r}$ is continuous in $\mathbb{C}$.

###### Proof.

Let us fix $r\in\mathbb{N}$ and consider any $z,h\in\mathbb{C}$. Then

$\displaystyle\bigl{|}F_{r}(z+h)-F_{r}(z)\bigr{|}$
$\displaystyle\leq\mathbb{E}\Bigl{|}|f_{r,z+h}(\tau)|^{2}-|f_{r,z}(\tau)|^{2}\Bigl{|}$
$\displaystyle\leq\mathbb{E}\Bigl{(}|f_{r,z+h}(\tau)-f_{r,z}(\tau)|\cdot\bigr{(}|f_{r,z+h}(\tau)|+|f_{r,z}(\tau)|\bigl{)}\Bigr{)}.$
(23)

We note that

$\displaystyle|f_{r,z+h}(\tau)|$
$\displaystyle\leq\sum_{j=0}^{r+1}\frac{|z+h|^{j}}{j!}+\frac{|z+h|^{r+2}}{(r+1)!}\leq\sum_{j=0}^{r+1}\frac{(|z|+|h|)^{j}}{j!}+\frac{(|z|+|h|)^{r+2}}{(r+1)!}=:\alpha_{z}(h)$

with probability $1$. Inserting this bound into (23) yields

$\displaystyle\bigl{|}F_{r}(z+h)-$ $\displaystyle
F_{r}(z)\bigr{|}\leq\bigl{(}\alpha_{z}(h)+\alpha_{z}(0)\bigr{)}\cdot\mathbb{E}|f_{r,z+h}(\tau)-f_{r,z}(\tau)|$
$\displaystyle=|h|\cdot\bigl{(}\alpha_{z}(h)+\alpha_{z}(0)\bigr{)}\cdot\mathbb{E}\Bigl{|}\sum_{j=1}^{r+1}\sum_{k=0}^{j-1}\frac{(z+h)^{j-1-k}z^{k}}{j!}+\tau^{r+1}\sum_{k=0}^{r+1}\frac{(z+h)^{r+1-k}z^{k}}{(r+1)!}\Bigr{|}$
$\displaystyle\leq|h|\cdot\bigl{(}\alpha_{z}(h)+\alpha_{z}(0)\bigr{)}\cdot\Bigl{(}\sum_{j=1}^{r+1}\frac{(|z|+|h|)^{j-1}}{(j-1)!}+(r+2)\cdot\frac{(|z|+|h|)^{r+1}}{(r+1)!}\Bigr{)}.$

The last expression tends to $0$ when $h\to 0$, which completes the proof. ∎

######  Lemma 3.

Let $k\in\mathbb{Z}_{+}$ and $\alpha\in[0,1]$. Then for each $\varepsilon>0$
there exists $\delta>0$ such that for any $z\in\mathbb{C}$ with
$|z-\alpha|<\delta$ we have

$\int\limits_{A_{\delta}}\bigl{|}\ln|t^{k}-z|\bigr{|}\,\mathrm{d}t<\varepsilon,$

where $A_{\delta}=\\{t\in[0,1]\colon|t^{k}-\alpha|\leq\delta\\}$.

###### Proof.

Let $k\in\mathbb{Z}_{+}$, $\alpha\in[0,1]$ and $\varepsilon>0$ be fixed. Let
us take into consideration only $\delta\in(0,\frac{1}{2})$. Then
$|t^{k}-z|\leq|t^{k}-\alpha|+|z-\alpha|<2\delta<1$ for all $z\in\mathbb{C}$
such that $|z-\alpha|<\delta$ and all $t\in A_{\delta}$. As a result,

$\int\limits_{A_{\delta}}\bigl{|}\ln|t^{k}-z|\bigr{|}\,\mathrm{d}t=-\int\limits_{A_{\delta}}\ln|t^{k}-z|\,\mathrm{d}t.$
(24)

Moreover,

$0\geq\int\limits_{A_{\delta}}\ln|t^{k}-z|\,\mathrm{d}t\geq\int\limits_{A_{\delta}}\ln|t^{k}-\Re(z)|\,\mathrm{d}t\geq\inf_{\begin{subarray}{c}x\in\mathbb{R}\\\
|x-\alpha|<\delta\end{subarray}}\ \
\int\limits_{A_{\delta}}\ln|t^{k}-x|\,\mathrm{d}t$ (25)

because
$|t^{k}-z|=\sqrt{\Re^{2}(t^{k}-z)+\Im^{2}(t^{k}-z)}\geq|\Re(t^{k}-z)|=|t^{k}-\Re(z)|$
and similarly $|\alpha-\Re(z)|\leq|\alpha-z|<\delta$. For $x\leq 0$ we obtain

$\int\limits_{A_{\delta}}\ln|t^{k}-x|\,\mathrm{d}t\geq
k\int\limits_{A_{\delta}}\ln t\,\mathrm{d}t>-\varepsilon$ (26)

if $\delta$ is sufficiently close to $0$. Now let us consider
$x\in(0,\alpha+\delta)$ and define $y=\sqrt[k]{x}$,

$d(\delta)=\left\\{\begin{array}[]{ll}0,&\text{when }\alpha-\delta\leq 0,\\\
\sqrt[k]{\alpha-\delta},&\text{when }\alpha-\delta>0.\end{array}\right.$

Then

$\displaystyle\int\limits_{A_{\delta}}\ln|t^{k}-x|\,\mathrm{d}t$
$\displaystyle=\int\limits_{d(\delta)}^{\sqrt[k]{\alpha+\delta}}\ln|t-y|\,\mathrm{d}t+\int\limits_{d(\delta)}^{\sqrt[k]{\alpha+\delta}}\ln\Bigl{(}\sum_{j=0}^{k-1}t^{k-1-j}y^{j}\Bigr{)}\,\mathrm{d}t$
$\displaystyle\geq\int\limits_{d(\delta)-y}^{\sqrt[k]{\alpha+\delta}-y}\ln|t|\,\mathrm{d}t+(k-1)\int\limits_{d(\delta)}^{\sqrt[k]{\alpha+\delta}}\ln
t\,\mathrm{d}t$ $\displaystyle\geq\ \
k\\!\\!\\!\\!\int\limits_{-\frac{1}{2}(\sqrt[k]{\alpha+\delta}-d(\delta))}^{\frac{1}{2}(\sqrt[k]{\alpha+\delta}-d(\delta))}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\ln|t|\,\mathrm{d}t>-\varepsilon$
(27)

for $\delta$ sufficiently close to $0$. In the last two lines of (27) we used
the fact that

$\int\limits_{a}^{b}\ln|t|\,\mathrm{d}t\geq\int\limits_{-\frac{b-a}{2}}^{\frac{b-a}{2}}\ln|t|\,\mathrm{d}t$

for all $a,b\in\mathbb{R}$ such that $a<b$; moreover, the last integral tends
to $0$ when $\frac{b-a}{2}\to 0$. By (24), (25), (26) and (27) we get the
desired claim. ∎

The following Lemma 4 is a generalization of Proposition 1 from [1].

######  Lemma 4.

For each $r\in\mathbb{N}$, the function $G_{r}$ is continuous in $\mathbb{C}$.

###### Proof of Lemma 4.

To see that $G_{r}$ is continuous in $0$, let us observe that

$0<1-\sum_{j=1}^{r+1}\frac{|z|^{j}}{j!}-\frac{|z|^{r+2}}{(r+1)!}\leq\bigl{|}f_{r,z}(t)\bigr{|}\leq
1+\sum_{j=1}^{r+1}\frac{|z|^{j}}{j!}+\frac{|z|^{r+2}}{(r+1)!}$

for $|z|$ sufficiently close to $0$ and for all $t\in[0,1]$. Hence,

$\ln\Bigl{(}1-\sum_{j=1}^{r+1}\frac{|z|^{j}}{j!}-\frac{|z|^{r+2}}{(r+1)!}\Bigr{)}\leq
G_{r}(z)\leq\ln\Bigl{(}1+\sum_{j=1}^{r+1}\frac{|z|^{j}}{j!}+\frac{|z|^{r+2}}{(r+1)!}\Bigr{)},$

which implies that $\displaystyle\lim_{z\to 0}G_{r}(z)=0=G_{r}(0)$.

Let us consider $z\in\mathbb{C}\setminus\\{0\\}$. Then $G_{r}(z)$ can be
expressed as

$G_{r}(z)=\ln\Bigl{(}\frac{|z|^{r+2}}{(r+1)!}\Bigr{)}+H_{r}(h_{r}(z)),$ (28)

where

$\displaystyle h_{r}\colon$ $\displaystyle\mathbb{C}\setminus\\{0\\}\ni
z\mapsto-\sum_{j=0}^{r+1}\frac{(r+1)!}{j!\cdot z^{r+2-j}}\in\mathbb{C},$ (29)
$\displaystyle H_{r}\colon$ $\displaystyle\mathbb{C}\ni
z\mapsto\int\limits_{0}^{1}\ln|t^{r+1}-z|\,\mathrm{d}t\in\mathbb{R}.$ (30)

To complete the proof, it suffices to show that $H_{r}$ is continuous in
$\mathbb{C}$.

Firstly, let us consider a fixed $z\in\mathbb{C}\setminus[0,1]$. Let us define

$\delta=\min_{t\in[0,1]}|z-t^{r+1}|$

and take any $(z_{n})_{n=1}^{\infty}\subset\mathbb{C}$ such that $z_{n}\to z$,
$z_{n}\neq z$ and $|z_{n}-z|<\frac{\delta}{2}$ for $n\in\mathbb{Z}_{+}$. Then
for all $t\in[0,1]$ we have
$|z_{n}-t^{r+1}|\geq|z-t^{r+1}|-|z_{n}-z|>\frac{\delta}{2}$. By Lagrange’s
mean value theorem and the triangle inequality, we obtain

$\left|\frac{\ln|t^{r+1}-z|-\ln|t^{r+1}-z_{n}|}{z-z_{n}}\right|\leq\left|\frac{\ln|t^{r+1}-z|-\ln|t^{r+1}-z_{n}|}{|t^{r+1}-z|-|t^{r+1}-z_{n}|}\right|=\frac{1}{\xi(t,n)}$

for some $\xi(t,n)$ falling between $|t^{r+1}-z|$ and $|t^{r+1}-z_{n}|$,
provided that $|t^{r+1}-z|\neq|t^{r+1}-z_{n}|$. Since both these numbers are
greater than $\frac{\delta}{2}$, we have $\xi(t,n)>\frac{\delta}{2}$ and thus

$\bigl{|}\ln|t^{r+1}-z|-\ln|t^{r+1}-z_{n}|\bigr{|}<\frac{2}{\delta}\cdot|z-z_{n}|$
(31)

for all $t\in[0,1]$. Note that the above inequality holds also when
$|t^{r+1}-z|=|t^{r+1}-z_{n}|$. From (31) it follows that
$t\mapsto\ln|t^{r+1}-z_{n}|$ converges uniformly to $t\mapsto\ln|t^{r+1}-z|$
for $t\in[0,1]$. Hence,

$\lim_{n\to\infty}\int\limits_{0}^{1}\ln|t^{r+1}-z_{n}|\,\mathrm{d}t=\int\limits_{0}^{1}\ln|t^{r+1}-z|\,\mathrm{d}t$

and continuity of $G_{r}$ in each point $z\in\mathbb{C}\setminus[0,1]$ is
proven.

Now let us consider a fixed $z\in[0,1]$ and set $\varepsilon>0$. By Lemma 3,
there exists $\delta>0$ such that for any $\zeta\in\mathbb{C}$ with
$|\zeta-z|<\delta$ we have

$\int\limits_{A_{\delta}}\bigl{|}\ln|t^{r+1}-\zeta|\bigr{|}\,\mathrm{d}t<\frac{\varepsilon}{3},$
(32)

where $A_{\delta}=\\{t\in[0,1]\colon|t^{r+1}-z|\leq\delta\\}$. Let us consider
any $(z_{n})_{n=1}^{\infty}\subset\mathbb{C}$ such that $z_{n}\to z$,
$z_{n}\neq z$ and $|z_{n}-z|<\frac{\delta}{2}$ for all $n\in\mathbb{Z}_{+}$.
For each $t\in[0,1]\setminus A_{\delta}$ we have $|t^{r+1}-z|>\delta$ and
$|t^{r+1}-z_{n}|>\frac{\delta}{2}$, $n\in\mathbb{Z}_{+}$. Thus, we can prove
the uniform convergence of $t\mapsto\ln|t^{r+1}-z_{n}|$ to
$t\mapsto\ln|t^{r+1}-z|$ for $t\in[0,1]\setminus A_{\delta}$ in a similar
fashion as in the case of $z\in\mathbb{C}\setminus[0,1]$, see (31). Hence, for
sufficiently big $n$ we obtain

$I=\Bigl{|}\ \int\limits_{[0,1]\setminus
A_{\delta}}\ln|t^{r+1}-z_{n}|\,\mathrm{d}t-\int\limits_{[0,1]\setminus
A_{\delta}}\ln|t^{r+1}-z|\,\mathrm{d}t\Bigr{|}<\frac{\varepsilon}{3}.$ (33)

By (32) and (33),

$\Bigl{|}\int\limits_{0}^{1}\ln|t^{r+1}-z_{n}|\,\mathrm{d}t-\int\limits_{0}^{1}\ln|t^{r+1}-z|\,\mathrm{d}t\Bigr{|}\leq
I+\Bigl{|}\int\limits_{A}\ln|t^{r+1}-z_{n}|\,\mathrm{d}t\Bigr{|}+\Bigl{|}\int\limits_{A}\ln|t^{r+1}-z|\,\mathrm{d}t\Bigr{|}<\varepsilon$

for sufficiently big $n$. This concludes the proof. ∎

## References

  * [1] T. Bochacik, M. Goćwin, P. M. Morkisz, and P. Przybyłowicz, _Randomized Runge-Kutta method – Stability and convergence under inexact information_ , J. Complex. 65 (2021), 101554 (21 pages).
  * [2] T. Bochacik and P. Przybyłowicz, _On the randomized Euler schemes for ODEs under inexact information_ , Numer. Algorithms (2022), https://doi.org/10.1007/s11075-022-01299-7.
  * [3] T. Daun, _On the randomized solution of initial value problems_ , J. Complex. 27 (2011), 300–311.
  * [4] T. Daun and S. Heinrich, _Complexity of parametric initial value problems in Banach spaces_ , J. Complex. 30 (2014), 392–429.
  * [5] E. Hairer and G. Wanner, _Solving Ordinary Differential Equations II, Stiff and Differential- Algebraic Problems, 2nd ed._ , Springer-Verlag, Berlin, 1996.
  * [6] M. Eisenmann, M. Kovács, R. Kruse, and S. Larsson, _On a randomized backward Euler method for nonlinear evolution equations with time-irregular coefficients_ , Found. Comp. Math. 19 (2019), 1387–1430.
  * [7] M. Goćwin and M. Szczęsny, _R_ andomized and quantum algorithms for solving initial-value problems in ordinary differential equations of order $k$, Opusc. Math. 28 (2008), 247–277.
  * [8] S. Heinrich and B. Milla, _The randomized complexity of initial value problems_ , J. Complex. 24 (2008), 77–88.
  * [9] D.J. Higham, Mean-square and asymptotic stability of the stochastic theta method, Siam J. Numer. Anal. 38 (2000), 753–769.
  * [10] A. Jentzen and A. Neuenkirch, _A random Euler scheme for Carathéodory differential equations_ , J. Comput. Appl. Math. 224 (2009), 346–359.
  * [11] B. Kacewicz, _Almost optimal solution of initial-value problems by randomized and quantum algorithms_ , J. Complex. 22 (2006), 676–690.
  * [12] B. Kacewicz, _Improved bounds on the randomized and quantum complexity of initial-value problems_ , J. Complex. 21 (2005), 740–756.
  * [13] B. Kacewicz, _Randomized and quantum algorithms yield a speed-up for initial-value problems_ , J. Complex. 20 (2004) 821–834.
  * [14] R. Kruse and Y. Wu, _Error Analysis of Randomized Runge–Kutta Methods for Differential Equations with Time-Irregular Coefficients_ , Comput. Methods Appl. Math. 17 (2017), 479–498.
  * [15] T. Mitsui and Y. Saito, _Stability Analysis of Numerical Schemes for Stochastic Differential Equations_ , SIAM J. Numer. Anal., 33 (1996), 2254–2267.
  * [16] G. Stengle, _Numerical methods for systems with measurable coefficients_ , Appl. Math. Lett. 3 (1990), 25–29.
  * [17] G. Stengle, _Error analysis of a randomized numerical method_ , Numer. Math. 70 (1995), 119–128.


# Thermodynamics of interacting holographic dark energy
with apparent horizon as an IR cutoff

Ahmad Sheykhi 111<EMAIL_ADDRESS>Department of Physics, Shahid Bahonar
University, P.O. Box 76175, Kerman, Iran
Research Institute for Astronomy and Astrophysics of Maragha (RIAAM), Maragha,
Iran

###### Abstract

As soon as an interaction between holographic dark energy and dark matter is
taken into account, the identification of IR cutoff with Hubble radius
$H^{-1}$, in flat universe, can simultaneously drive accelerated expansion and
solve the coincidence problem. Based on this, we demonstrate that in a non-
flat universe the natural choice for IR cutoff could be the apparent horizon
radius, $\tilde{r}_{A}={1}/{\sqrt{H^{2}+k/a^{2}}}$. We show that any
interaction of dark matter with holographic dark energy, whose infrared cutoff
is set by the apparent horizon radius, implies an accelerated expansion and a
constant ratio of the energy densities of both components thus solving the
coincidence problem. We also verify that for a universe filled with dark
energy and dark matter the Friedmann equation can be written in the form of
the modified first law of thermodynamics, $dE=T_{h}dS_{h}+WdV$, at apparent
horizon. In addition, the generalized second law of thermodynamics is
fulfilled in a region enclosed by the apparent horizon. These results hold
regardless of the specific form of dark energy and interaction term. Our study
might reveal that in an accelerating universe with spatial curvature, the
apparent horizon is a physical boundary from the thermodynamical point of
view.

##  I Introduction

The combined analysis of cosmological observations reveal that nearly three
quarters of our universe consists of a mysterious energy component usually
dubbed “dark energy” which is responsible for the cosmic expansion, and the
remaining part consists of pressureless matter Rie . The nature of such
previously unforeseen energy still remains a complete mystery, except for the
fact that it has negative pressure. In this new conceptual set up, one of the
important questions concerns the thermodynamical behavior of the accelerated
expanding universe driven by dark energy. It is important to ask whether
thermodynamics in an accelerating universe can reveal some properties of dark
energy. The profound connection between thermodynamics and gravity has been
observed in the cosmological situations Cai2 ; Cai3 ; CaiKim ; Fro ; Wang ;
Cai4 ; Shey1 ; Shey2 ; Shey3 . This connection implies that the
thermodynamical properties can help understand the dark energy, which gives
strong motivation to study thermodynamics in the accelerating universe. It is
also of great interest to investigate the validity of the generalized second
law of thermodynamics in the accelerating universe driven by dark energy
wangb0 . The generalized second law of thermodynamics is an important
principle in governing the development of the nature.

An interesting attempt for probing the nature of dark energy within the
framework of quantum gravity, is the so-called “Holographic Dark Energy” (HDE)
proposal. This model which has arisen a lot of enthusiasm recently Coh ; Li ;
Huang ; Hsu ; HDE ; Setare1 ; wang0 ; wang1 , is motivated from the
holographic hypothesis Suss1  and has been tested and constrained by various
astronomical observations Xin . It is important to note that in the
literature, various scenarios of HDE have been studied via considering
different system’s IR cutoff. In the absence of interaction between dark
matter and dark energy in flat universe, Li Li  discussed three choices for
the length scale $L$ which is supposed to provide an IR cutoff. The first
choice is the Hubble radius, $L=H^{-1}$ Hsu , which leads to a wrong equation
of state, namely that for dust. The second option is the particle horizon
radius. In this case it is impossible to obtain an accelerated expansion. Only
the third choice, the identification of $L$ with the radius of the future
event horizon gives the desired result, namely a sufficiently negative
equation of state to obtain an accelerated universe.

However, as soon as an interaction between dark energy and dark matter is
taken into account, the first choice, $L=H^{-1}$, in flat universe, can
simultaneously drive accelerated expansion and solve the coincidence problem
pav1 . Based on this, we demonstrate that in a non-flat universe the natural
choice for IR cutoff could be the apparent horizon radius. We show that any
interaction of pressureless dark matter with HDE, whose infrared cutoff is set
by the apparent horizon radius, implies a constant ratio of the energy
densities of both components thus solving the coincidence problem. Besides, it
was argued that for an accelerating universe inside the event horizon the
generalized second law does not satisfy, while the accelerating universe
enveloped by the Hubble horizon satisfies the generalized second law Jia .
This implies that the event horizon in an accelerating universe might not be a
physical boundary from the thermodynamical point of view. Thus, it looks that
we need to define a convenient horizon that satisfies all of our accepted
principles in a universe with any spacial curvature. In the next section, we
study the interacting HDE with apparent horizon as an IR cutoff. In section
III, we examine the first law of thermodynamics on the apparent horizon in an
accelerating universe with spacial curvature. In section IV, we investigate
the validity of the generalized second law of thermodynamics in a region
enclosed by the apparent horizon. The last section is devoted to conclusions.

##  II Interacting HDE with apparent horizon as an IR cutoff

We consider a homogenous and isotropic Friedmann-Robertson-Walker (FRW)
universe which is described by the line element

$ds^{2}={h}_{\mu\nu}dx^{\mu}dx^{\nu}+\tilde{r}^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2}),$ (1)

where $\tilde{r}=a(t)r$, $x^{0}=t,x^{1}=r$, the two dimensional metric
$h_{\mu\nu}$=diag $(-1,a^{2}/(1-kr^{2}))$. Here $k$ denotes the curvature of
space with $k=0,1,-1$ corresponding to open, flat, and closed universes,
respectively. A closed universe with a small positive curvature
($\Omega_{k}\simeq 0.01$) is compatible with observations spe . Then, the
dynamical apparent horizon, a marginally trapped surface with vanishing
expansion, is determined by the relation
$h^{\mu\nu}\partial_{\mu}\tilde{r}\partial_{\nu}\tilde{r}=0$, which implies
that the vector $\nabla\tilde{r}$ is null on the apparent horizon surface. The
apparent horizon was argued as a causal horizon for a dynamical spacetime and
is associated with gravitational entropy and surface gravity Hay2 ; Bak . A
simple calculation gives the apparent horizon radius for the FRW universe

$\tilde{r}_{A}=\frac{1}{\sqrt{H^{2}+k/a^{2}}}.$ (2)

The corresponding Friedmann equation takes the form

$\displaystyle H^{2}+\frac{k}{a^{2}}=\frac{8\pi
G}{3}\left(\rho_{m}+\rho_{D}\right),$ (3)

where $\rho_{m}$ and $\rho_{D}$ are the energy density of dark matter and dark
energy inside apparent horizon, respectively. Since we consider the
interaction between dark matter and dark energy, $\rho_{m}$ and $\rho_{D}$ do
not conserve separately; they must rather enter the energy balances

$\displaystyle\dot{\rho}_{m}+3H\rho_{m}=Q,$ (4)
$\displaystyle\dot{\rho}_{D}+3H\rho_{D}(1+w_{D})=-Q.$ (5)

where $w_{D}=p_{D}/\rho_{D}$ is the equation of state parameter of HDE, and
$Q$ stands for the interaction term. We also ignore the baryonic matter
($\Omega_{BM}\approx 0.04$) in comparison with dark matter and dark energy
($\Omega_{DM}+\Omega_{DE}\approx 0.96$). We shall assume the ansatz
$Q=\Gamma\rho_{D}$ with $\Gamma>0$ which means that there is an energy
transfer from the dark energy to dark matter. It is important to note that the
continuity equations imply that the interaction term should be a function of a
quantity with units of inverse of time (a first and natural choice can be the
Hubble factor $H$) multiplied with the energy density. Therefore, the
interaction term could be in any of the following forms: (i) $Q\propto
H\rho_{D}$, (ii) $Q\propto H\rho_{m}$, or (iii) $Q\propto
H(\rho_{m}+\rho_{D})$. However, we can present the above three forms in one
expression as $Q=\Gamma\rho_{D}$, where

$\displaystyle\begin{array}[]{ll}\Gamma=3b^{2}H\hskip 36.98866pt{\rm for}\ \
Q\propto H\rho_{D},&\\\ \Gamma=3b^{2}Hu\hskip 31.2982pt{\rm for}\ \ Q\propto
H\rho_{m},&\\\ \Gamma=3b^{2}H(1+u)\ \ {\rm for}\ \ Q\propto
H(\rho_{m}+\rho_{D}),&\end{array}$ (9)

with $b^{2}$ is a coupling constant and $u=\rho_{m}/\rho_{D}$ is the ratio of
energy densities. The freedom of choosing the specific form of the interaction
term $Q$ stems from our incognizance of the origin and nature of dark energy
as well as dark matter. Moreover, a microphysical model describing the
interaction between the dark components of the universe is not available
nowadays. If we introduce, as usual, the fractional energy densities such as

$\displaystyle\Omega_{m}=\frac{8\pi G\rho_{m}}{3H^{2}},\hskip
14.22636pt\Omega_{D}=\frac{8\pi G\rho_{D}}{3H^{2}},\hskip
14.22636pt\Omega_{k}=\frac{k}{H^{2}a^{2}},$ (10)

then, the Friedmann equation can be written as
$\Omega_{m}+\Omega_{D}=1+\Omega_{k}$. In terms of the apparent horizon radius,
we can rewrite the Friedmann equation as

$\frac{1}{\tilde{r}_{A}^{2}}=\frac{8\pi G}{3}\left(\rho_{m}+\rho_{D}\right).$
(11)

For completeness, we give the deceleration parameter

$q=-\frac{\ddot{a}}{aH^{2}}=-1-\frac{\dot{H}}{H^{2}},$ (12)

which combined with the Hubble parameter and the dimensionless density
parameters form a set of useful parameters for the description of the
astrophysical observations. It is a matter of calculation to show that

$q=-(1+\Omega_{k})+\frac{3}{2}\Omega_{D}(1+u+w_{D}).$ (13)

The evolution of $u$ is governed by

$\dot{u}=3Hu\left[w_{D}+\frac{1+u}{u}\frac{\Gamma}{3H}\right].$ (14)

We assume the HDE density has the form

$\rho_{D}=\frac{3c^{2}}{8\pi G\tilde{r}_{A}^{2}},$ (15)

where $c^{2}$ is a constant, the coefficient $3$ is for convenient, and we
have set the apparent horizon radius $L={\tilde{r}_{A}}$ as system’s IR cutoff
in holographic model of dark energy. Inserting Eq. (15) in Eq. (11)
immediately yields

$\rho_{m}=\frac{3(1-c^{2})}{8\pi G\tilde{r}_{A}^{2}}.$ (16)

Thus we reach

$u=\frac{\rho_{m}}{\rho_{D}}=\frac{1-c^{2}}{c^{2}}.$ (17)

This implies that the ratio of the energy densities is a constant; thus the
coincidence problem can be solved. Taking the derivative of Eq. (15) we get

$\dot{\rho}_{D}=-2\rho_{D}\frac{\dot{\tilde{r}_{A}}}{\tilde{r}_{A}}=-3c^{2}H\rho_{D}(1+u+w_{D}).$
(18)

where we have employed Eqs. (4), (5) and (11). Combining this equation with
(5) we obtain

$w_{D}=-\left(1+\frac{1}{u}\right)\frac{\Gamma}{3H}.$ (19)

Substituting $w_{D}$ into (13), we find

$q=-(1+\Omega_{k})-\frac{3}{2}\Omega_{D}(1+u)\left(\frac{\Gamma}{3Hu}-1\right).$
(20)

The interaction parameter $\frac{\Gamma}{3H}$ together with the energy density
ratio $u$ determine the equation of state parameter. In the absence of
interaction, we encounter dust with $w_{D}=0$. For the choice
$L=\tilde{r}_{A}$ an interaction is the only way to have an equation of state
different from that for dust. Any decay of the dark energy component into
pressureless matter is necessarily accompanied by an equation of state
$w_{D}<0$. The existence of an interaction has another interesting
consequence. Inserting expression $w_{D}$ into (14) leads to $\dot{u}=0$,
i.e., $u$ = const. Thus, any interaction of dark matter with HDE, whose
infrared cutoff is set by the apparent horizon radius, implies an accelerated
expansion and a constant ratio of the energy densities, irrespective of the
specific structure of the interaction. It is important to note that although
choosing $L=H^{-1}$, in a spatially flat universe, can drive accelerated
expansion and solve the coincidence problem pav1 , but taking into account the
spatial curvature term gives rise to an additional dynamics which implies a
small (compared with the Hubble rate) change of the energy density ratio; thus
the coincidence problem cannot be solved exactly (see pav2  for details). This
implies that in an accelerating universe with spacial curvature the Hubble
radius $H^{-1}$ is not a convenient choice.

In summary, in a universe with spacial curvature, the identification of IR
cutoff with apparent horizon radius $\tilde{r}_{A}$ is not only the most
obvious but also the simplest choice which can simultaneously drive
accelerated expansion and solve the coincidence problem. It is important to
note that the interaction is essential to simultaneously solve the coincidence
problem and have late acceleration. There is no non-interacting limit, since
in the absence of interaction, i.e., $\Gamma=0$, there is no acceleration.

##  III First law of thermodynamics

In this section we are going to examine the first law of thermodynamics. In
particular, we show that for a closed universe filled with HDE and dark matter
the Friedmann equation can be written directly in the form of the modified
first law of thermodynamics at apparent horizon regardless of the specific
form of the dark energy. The associated temperature with the apparent horizon
can be defined as $T=\kappa/2\pi$, where $\kappa$ is the surface gravity
$\kappa=\frac{1}{\sqrt{-h}}\partial_{\mu}\left(\sqrt{-h}h^{\mu\nu}\partial_{\mu\nu}\tilde{r}\right).$
Then one can easily show that the surface gravity at the apparent horizon of
FRW universe can be written as

$\kappa=-\frac{1}{\tilde{r}_{A}}\left(1-\frac{\dot{\tilde{r}}_{A}}{2H\tilde{r}_{A}}\right).$
(21)

When $\dot{\tilde{r}}_{A}\leq 2H\tilde{r}_{A}$, the surface gravity
$\kappa\leq 0$, which leads the temperature $T\leq 0$ if one defines the
temperature of the apparent horizon as $T=\kappa/2\pi$ . Physically it is not
easy to accept the negative temperature, the temperature on the apparent
horizon should be defined as $T=|\kappa|/2\pi$. Recently the connection
between temperature on the apparent horizon and the Hawking radiation has been
considered in cao , which gives more solid physical implication of the
temperature associated with the apparent horizon.

Taking differential form of equation (11) and using Eqs. (4) and (5), we can
get the differential form of the Friedmann equation

$\frac{1}{4\pi
G}\frac{d\tilde{r}_{A}}{\tilde{r}_{A}^{3}}=H\rho_{D}\left(1+u+w_{D}\right)dt.$
(22)

Multiplying both sides of the equation (22) by a factor
$4\pi\tilde{r}_{A}^{3}\left(1-\frac{\dot{\tilde{r}}_{A}}{2H\tilde{r}_{A}}\right)$,
and using the expression (21) for the surface gravity, after some
simplification one can rewrite this equation in the form

$\displaystyle-\frac{\kappa}{2\pi}\frac{2\pi\tilde{r}_{A}d\tilde{r}_{A}}{G}$
$\displaystyle=$ $\displaystyle
4\pi\tilde{r}_{A}^{3}H\rho_{D}\left(1+u+w_{D}\right)dt-2\pi\tilde{r}_{A}^{2}\rho_{D}\left(1+u+w_{D}\right)d\tilde{r}_{A}.$
(23)

$E=(\rho_{m}+\rho_{D})V$ is the total energy content of the universe inside a
$3$-sphere of radius $\tilde{r}_{A}$, where
$V=\frac{4\pi}{3}\tilde{r}_{A}^{3}$ is the volume enveloped by 3-dimensional
sphere with the area of apparent horizon $A=4\pi\tilde{r}_{A}^{2}$. Taking
differential form of the relation
$E=(\rho_{m}+\rho_{D})\frac{4\pi}{3}\tilde{r}_{A}^{3}$ for the total matter
and energy inside the apparent horizon, we get

$dE=4\pi\tilde{r}_{A}^{2}(\rho_{m}+\rho_{D})d\tilde{r}_{A}+\frac{4\pi}{3}\tilde{r}_{A}^{3}(\dot{\rho}_{m}+\dot{\rho}_{D})dt.$
(24)

Using Eqs. (4) and (5), we obtain

$dE=4\pi\tilde{r}_{A}^{2}\rho_{D}(1+u)d\tilde{r}_{A}-4\pi\tilde{r}_{A}^{3}H\rho_{D}\left(1+u+w_{D}\right)dt.$
(25)

Substituting this relation into (23), and using the relation between
temperature and the surface gravity, we get the modified first law of
thermodynamics on the apparent horizon

$dE=T_{h}dS_{h}+WdV,$ (26)

where $S_{h}={A}/{4G}$ is the entropy associated to the apparent horizon, and

$W=\frac{1}{2}(\rho_{m}+\rho_{D}-p_{D})=\frac{1}{2}\rho_{D}\left(1+u-w_{D}\right)$
(27)

is the matter work density Hay2 . The work density term is regarded as the
work done by the change of the apparent horizon, which is used to replace the
negative pressure if compared with the standard first law of thermodynamics,
$dE=TdS-pdV$. For a pure de Sitter space, $\rho_{m}+\rho_{D}=-p_{D}$, then our
work term reduces to the standard $-p_{D}dV$ and we obtain exactly the first
law of thermodynamics.

##  IV Generalized Second law of thermodynamics

In this section we turn to investigate the validity of the generalized second
law of thermodynamics in a region enclosed by the apparent horizon.
Differentiating Eq. (11) with respect to the cosmic time and using Eqs. (4)
and (5) we get

$\dot{\tilde{r}}_{A}=4\pi GH{\tilde{r}_{A}^{3}}\rho_{D}(1+u+w_{D}).$ (28)

One can see from the above equation that $\dot{\tilde{r}}_{A}>0$ provided
condition $w_{D}>-1-u$, holds. Let us now turn to find out $T_{h}\dot{S_{h}}$:

$T_{h}\dot{S_{h}}=\frac{1}{2\pi\tilde{r}_{A}}\left(1-\frac{\dot{\tilde{r}}_{A}}{2H\tilde{r}_{A}}\right)\frac{d}{dt}\left(\frac{\pi\tilde{r}_{A}^{2}}{G}\right).$
(29)

After some simplification and using Eq. (28) we get

$T_{h}\dot{S_{h}}=4\pi
H{\tilde{r}_{A}^{3}}\rho_{D}(1+u+w_{D})\left(1-\frac{\dot{\tilde{r}}_{A}}{2H\tilde{r}_{A}}\right).$
(30)

As we argued above the term
$\left(1-\frac{\dot{\tilde{r}}_{A}}{2H\tilde{r}_{A}}\right)$ is positive to
ensure $T_{h}>0$, however, in an accelerating universe the equation of state
parameter of dark energy may cross the phantom divide, i.e., $w_{D}<-1-u$.
This indicates that the second law of thermodynamics, $\dot{S_{h}}\geq 0$,
does not hold on the apparent horizon. Then the question arises, “will the
generalized second law of thermodynamics,
$\dot{S_{h}}+\dot{S_{m}}+\dot{S_{D}}\geq 0$, can be satisfied in a region
enclosed by the apparent horizon?” The entropy of dark energy plus dark matter
inside the apparent horizon, $S=S_{m}+S_{D}$, can be related to the total
energy $E=(\rho_{m}+\rho_{D})V$ and pressure $p_{D}$ in the horizon by the
Gibbs equation Pavon2

$TdS=d[(\rho_{m}+\rho_{D})V]+p_{D}dV=V(d\rho_{m}+d\rho_{D})+\rho_{D}(1+u+w_{D})dV,$
(31)

where $T=T_{m}=T_{D}$ and $S=S_{m}+S_{D}$ are the temperature and the total
entropy of the energy and matter content inside the horizon, respectively.
Here we assumed that the temperature of both dark components are equal, due to
their mutual interaction. We also limit ourselves to the assumption that the
thermal system bounded by the apparent horizon remains in equilibrium so that
the temperature of the system must be uniform and the same as the temperature
of its boundary. This requires that the temperature $T$ of the energy content
inside the apparent horizon should be in equilibrium with the temperature
$T_{h}$ associated with the apparent horizon, so we have $T=T_{h}$Pavon2 .
This expression holds in the local equilibrium hypothesis. If the temperature
of the fluid differs much from that of the horizon, there will be spontaneous
heat flow between the horizon and the fluid and the local equilibrium
hypothesis will no longer hold. This is also at variance with the FRW
geometry. In general, when we consider the thermal equilibrium state of the
universe, the temperature of the universe is associated with the apparent
horizon. Therefore from the Gibbs equation (31) we can obtain

$T_{h}(\dot{S_{m}}+\dot{S_{D}})=4\pi{\tilde{r}_{A}^{2}}\rho_{D}(1+u+w_{D})\dot{\tilde{r}}_{A}-4\pi
H{\tilde{r}_{A}^{3}}\rho_{D}(1+u+w_{D}).$ (32)

To check the generalized second law of thermodynamics, we have to examine the
evolution of the total entropy $S_{h}+S_{m}+S_{D}$. Adding equations (30) and
(32), we get

$T_{h}(\dot{S}_{h}+\dot{S}_{m}+\dot{S}_{D})=2\pi{\tilde{r}_{A}^{2}}\rho_{D}(1+u+w_{D})\dot{\tilde{r}}_{A}=\frac{A}{2}\rho_{D}(1+u+w_{D})\dot{\tilde{r}}_{A}.$
(33)

where $A>0$ is the area of apparent horizon. Substituting
$\dot{\tilde{r}}_{A}$ from Eq. (28) into (33) we get

$T_{h}(\dot{S}_{h}+\dot{S}_{m}+\dot{S}_{D})=2\pi
GAH{\tilde{r}_{A}}^{3}\rho^{2}_{D}(1+u+w_{D})^{2}.$ (34)

The right hand side of the above equation cannot be negative throughout the
history of the universe, which means that
$\dot{S_{h}}+\dot{S_{m}}+\dot{S_{D}}\geq 0$ always holds. This indicates that
for a universe with spacial curvature filled with interacting dark components,
the generalized second law of thermodynamics is fulfilled in a region enclosed
by the apparent horizon.

##  V conculusions

It is worthwhile to note that in the literature, various scenarios of HDE have
been studied via considering different system’s IR cutoff. In the absence of
interaction the convenient choice for the IR cutoff are the radial size of the
horizon $R_{h}$ and the radius of the event horizon measured on the sphere of
the horizon $L=ar(t)$ in spatially flat and curved universe, respectively.
Although, in these cases the HDE gives the observation value of dark energy in
the universe and can drive the universe to an accelerated expansion phase, but
an obvious drawback concerning causality appears. Event horizon is a global
concept of spacetime; existence of event horizon of the universe depends on
future evolution of the universe; and event horizon exists only for universe
with forever accelerated expansion. However, as soon as an interaction between
dark energy and dark matter is taken into account, the identification of $L$
with $H^{-1}$ in flat universe, can simultaneously drive accelerated expansion
and solve the coincidence problem pav1 . The Hubble radius is not only the
most obvious but also the simplest choice in flat universe.

In this paper, we demonstrated that in a universe with spacial curvature the
natural choice for IR cutoff could be the apparent horizon radius,
$\tilde{r}_{A}={1}/{\sqrt{H^{2}+k/a^{2}}}$. We showed that any interaction of
pressureless dark matter with HDE, whose infrared cutoff is set by the
apparent horizon radius, implies a constant ratio of the energy densities of
both dark components thus solving the coincidence problem. In addition, we
examined the validity of the first and the generalized second law of
thermodynamics for a universe filled with mutual interacting dark components
in a region enclosed by the apparent horizon. These results hold regardless of
the specific form of dark energy and interaction term $Q$. Our study further
supports that in a universe with spatial curvature, the apparent horizon is a
physical boundary from the thermodynamical point of view.

###### Acknowledgements.

I thank the anonymous referee for constructive comments. I am also grateful to
Prof. B. Wang for helpful discussions and reading the manuscript. This work
has been supported by Research Institute for Astronomy and Astrophysics of
Maragha.

## References

  * (1) A.G. Riess, et al., Astron. J. 116 (1998) 1009;
S. Perlmutter, et al., Astrophys. J. 517 (1999) 565;
S. Perlmutter, et al., Astrophys. J. 598 (2003) 102;
P. de Bernardis, et al., Nature 404 (2000) 955.

  * (2) M. Akbar and R. G. Cai, Phys. Rev. D 75, 084003 (2007).
  * (3) R. G. Cai and L. M. Cao, Phys.Rev. D 75, 064008 (2007).
  * (4) R. G. Cai and S. P. Kim, JHEP 0502, 050 (2005).
  * (5) A. V. Frolov and L. Kofman, JCAP 0305, 009 (2003);
U. K. Danielsson, Phys. Rev. D 71, 023516(2005) ;
R. Bousso, Phys. Rev. D 71, 064024 (2005);
G. Calcagni, JHEP 0509, 060 (2005).

  * (6) B. Wang, E. Abdalla and R. K. Su, Phys.Lett. B 503, 394 (2001);
B. Wang, E. Abdalla and R. K. Su, Mod. Phys. Lett. A 17, 23 (2002);
R. G. Cai and Y. S. Myung, Phys. Rev. D 67, 124021 (2003).

  * (7) R. G. Cai and L. M. Cao, Nucl. Phys. B 785 (2007) 135.
  * (8) A. Sheykhi, B. Wang and R. G. Cai, Nucl. Phys. B 779 (2007)1.
  * (9) A. Sheykhi, B. Wang and R. G. Cai, Phys. Rev. D 76 (2007) 023515.
  * (10) A. Sheykhi, B. Wang, Phys. Lett. B 678 (2009) 434.
  * (11) B. Wang, Y. Gong, E. Abdalla, Phys. Rev. D 74 (2006) 083520.
  * (12) A. Cohen, D. Kaplan, A. Nelson, Phys. Rev. Lett. 82 (1999) 4971\.
  * (13) M. Li, Phys. Lett. B 603 (2004) 1.
  * (14) Q. G. Huang, M. Li, JCAP 0408 (2004) 013.
  * (15) S. D. H. Hsu, Phys. Lett. B 594 (2004) 13.
  * (16) E. Elizalde, S. Nojiri, S.D. Odintsov, P. Wang, Phys. Rev. D 71 (2005) 103504;
B. Guberina, R. Horvat, H. Stefancic, JCAP 0505 (2005) 001;
B. Guberina, R. Horvat, H. Nikolic, Phys. Lett. B 636 (2006) 80;
H. Li, Z. K. Guo, Y. Z. Zhang, Int. J. Mod. Phys. D 15 (2006) 869;
Q. G. Huang, Y. Gong, JCAP 0408 (2004) 006;
J. P. B. Almeida, J. G. Pereira, Phys. Lett. B 636 (2006) 75;
Y. Gong, Phys. Rev. D 70 (2004) 064029;
B. Wang, E. Abdalla, R. K. Su, Phys. Lett. B 611 (2005) 21.

  * (17) M. R. Setare, S. Shafei, JCAP 09 (2006) 011;
M. R. Setare, Phys. Lett. B 644 (2007) 99;
M. R. Setare, JCAP 0701 (2007) 023;
M. R. Setare, Phys. Lett. B 654 (2007) 1;
M. R. Setare, Phys. Lett. B 642 (2006) 421.

  * (18) B. Wang, Y. Gong and E. Abdalla, Phys. Lett. B 624 (2005) 141;
B. Wang, C. Y. Lin. D. Pavon and E. Abdalla, Phys. Lett. B 662 (2008) 1;
A. Sheykhi, Phys. Lett. B 681 (2009) 205\.

  * (19) B. Wang, C. Y. Lin and E. Abdalla, Phys. Lett. B 637 (2005) 357;
M. R. Setare, Phys. Lett. B 642 (2006)1.

  * (20) G. t Hooft, gr-qc/9310026;
L. Susskind, J. Math. Phys. 36 (1995) 6377\.

  * (21) X. Zhang, F. Q. Wu, Phys. Rev. D 72 (2005) 043524;
X. Zhang, F. Q. Wu, Phys. Rev. D 76 (2007) 023502;
Q. G. Huang, Y.G. Gong, JCAP 0408 (2004) 006;
K. Enqvist, S. Hannestad, M. S. Sloth, JCAP 0502 (2005) 004;
J. Y. Shen, B. Wang, E. Abdalla, R.K. Su, Phys. Lett. B 609 (2005) 200.

  * (22) D. Pavon, W. Zimdahl, Phys. Lett. B 628 (2005) 206\.
  * (23) J. Zhou, B. Wang, Y. Gong, E. Abdalla, Phys. Lett. B 652 (2007) 86.
  * (24) D. N. Spergel, Astrophys. J. Suppl. 148 (2003) 175;
C. L. Bennett, et al., Astrophys. J. Suppl. 148 (2003) 1;
U. Seljak, A. Slosar, P. McDonald, JCAP 0610 (2006) 014;
D. N. Spergel, et al., Astrophys. J. Suppl. 170 (2007) 377.

  * (25) S.A. Hayward, S. Mukohyana, and M. C. Ashworth, Phys. Lett. A 256, 347 (1999);
S. A. Hayward, Class. Quantum Grav. 15, 3147 (1998).

  * (26) D. Bak and S. J. Rey, Class. Quantum Grav. 17, L83 (2000).
  * (27) L. Amendola, G. Camargo Campos, R. Rosenfeld, Phys. Rev. D 75 (2007) 083506;
Z.K. Guo, N. Ohta, S. Tsujikawa, Phys. Rev. D 76 (2007) 023508.

  * (28) W. Zimdahl, D. Pavon, Class. Quantum. Grav. 24 (2007) 5461.
  * (29) R.G. Cai, L.M. Cao, Y.P. Hu, Class. Quantum. Grav. 26 155018 (2009);
R. Li, J. R. Ren, D. F. Shi, Phys. Lett. B 670 (2009) 446.

  * (30) G. Izquierdo and D. Pavon, Phys.Lett. B 633 (2006) 420\.