# Invariant Measure for Neutral Stochastic Functional Differential Equations
with Non-Lipschitz Coeffiecients

Andriy Stanzhytskyi<EMAIL_ADDRESS>Department of Physics and
Mathematics, Igor Sikorsky Kyiv Polytechnic Institute, Ukraine  Oleksandr
Stanzhytskyi<EMAIL_ADDRESS>Department of Mathematics, Taras Shevchenko
National University of Kyiv, Ukraine  Oleksandr Misiats<EMAIL_ADDRESS>Department of Mathematics, Virginia Commonwealth University, Richmond, VA,
23284, USA

###### Abstract

In this work we study the long time behavior of nonlinear stochastic
functional-differential equations of neutral type in Hilbert spaces with non-
Lipschitz nonlinearities. We establish the existence of invariant measures in
the shift spaces for such equations. Our approach is based on Krylov-
Bogoliubov theorem on the tightness of the family of measures.

##  1 Introduction

In this work we study the asymptotic behaviour of the solutions of neutral
type stochastic functional-differential equations of the form

$\displaystyle d[u(t)+g(u_{t})]=[Au+f(u_{t})]dt+\sigma(u_{t})dW(t)\text{ for
}t>0;$ (1) $\displaystyle u(t)=\varphi(t),t\in[-h,0),\ h>0.$

Here $A$ is an inifinitesimal generator of a strong continuous semigroup
$\\{S(t),t\geq 0\\}$ of bounded linear operators in real separable Hilbert
space $H$. The noise $W(t)$ is a $Q$-Wiener process on a separable Hilbert
space $K$. For any $h>0$ denote $C_{h}:=C([-h,0],H)$ to be a space of
continuous $H$-valued functions $\varphi:[-h,0]\to H$, equipped with the norm

$\|\varphi\|_{C_{h}}:=\sup_{t\in[-h,0]}\|\varphi(t)\|_{H},$

where $\|\cdot\|_{H}$ stands for the norm in $H$. Throughout this paper,
$\|\cdot\|_{H}$ will be denoted with $\|\cdot\|$. The solution $u(t)$ of (1)
is sometimes referred as a state process. We also denote $u_{t}:=u(t+\theta)$,
where $\theta\in[-h,0)$. The functionals $f$ and $g$ map $C_{h}$ to $H$, and
$\sigma:C_{h}\to\mathcal{L}_{2}^{0}$, where
$\mathcal{L}_{2}^{0}=\mathcal{L}(Q^{1/2}K,H)$ is the space of Hilbert-Schmidt
operators from $Q^{1/2}K$ to $H$. Finally, $\varphi:[-h,0]\times\Omega\to H$
is the initial condition, where $(\Omega,\mathcal{F},P)$ is the probability
space.

Differential equations of neutral type arise in modeling the memory effects in
various dynamical processes (e.g. heat conduction). While the classic heat
equation $u_{t}=\Delta u$ accurately describes the heat transfer by conduction
in many materials, this model has two important drawbacks. First, the classic
heat equation does not take into account the memory effects which may be
present in some materials, especially at low temperatures. Most importantly,
the heat equation predicts an unrealistic result: that a thermal disturbance
at one point of the body is instantly felt everywhere in the body, although
not uniformly. In other words in Fourier heat conductors, finite thermal
discontinuities must propagate with infinite speed. These observations lead
one to believe that in materials with memory, the Fourier’s law of heat
conduction may be an approximation (perhaps for sufficiently steady
temperature fields) to a more general nonlinear constitutive assumption
relating the heat flux to the material’s thermal history. Gurtin and Pipkin
[10] have proposed a memory theory of heat conduction which is independent of
the present value of the temperature gradient, generalizes constitutive
relations deduced from kinetic theory by Maxwell [20] and Cattaneo [3], and
has finite heat propagation speeds. The linearized version of Gurkin-Pipkin’s
model reads as

$\dot{u}(x,t)+\int_{0}^{\infty}\beta(s)\dot{u}(x,t-s)\,ds=C\Delta u(x,t),$

which is a particular example of a neutral type differential equation. Similar
memory effects emerge in Hodgkin-Huxley model, Dawson-Fleming model of
population genetics [6] and other models.

The natural first step in the analysis of (1) is to establish the existence
and uniqueness of solutions and their continuous dependence on the initial
data. The question of well-posedness was addressed, e.g. in [18] under the
Lipschitz conditions. This condition was somewhat relaxed in [13]. The work
[1] addressed similar questions for the equations with impulsive action, as
well as with infinite delays. However, in our opinion, the proof of the main
result in aforementioned works is not quite complete. The reason is that this
proof is based on a successive approximation scheme, which is implicit. In
particular, the existence of a solution at each iteration step is highly
nontrivial and requires a separate proof.

On the other hand, [8] correctly establishes the existence and uniqueness of
mild solution under the Lipschits and linear growth assumptions. Under the
similar condition the works [26, 2] establish the existence and uniqueness of
solutions if the equation is perturbed with Brownian motion, as well as
fractional Brownian motion. However, the continuous dependence on the initial
conditions was not established there.

To the best of our knowledge, the question of the existence of invariant
measures for nonlinear equations of type (1) has not been addressed in
literature. One exception is the work [17], which established the existence of
a stationary solution for a linear equation of type (1).

The main goal of this paper is to establish the existence of invariant measure
for the equation (1) using Krylov-Bogoliubov theory on the tightness of a
family of measures [15]. More precisely, we will implement the compactness
approach of Da Prato and Zabczyk [5]. This approach was used in establishing
the existence of invariant measure for a large class of stochastic partial
differential equations with no delay [24, 23, 22, 12]. The results on long
time behavior of ordinary differential equations of similar nature were
obtained in [29, 28]. The existence of periodic solutions for neutral
functional differential equations in deterministic case were derived in [16].
Analogous results for non-neutral type stochastic functional differential
equations were established in [21].

In order to apply the compactness approach, we need to establish that the
solutions are both Markovian and Feller. These properties were established for
(1) in [14] in case when $A$ is Laplacian. The analog of comparison theorem
for such equations was established in [30].

The paper is organized as follows. In Section 2 we introduce the necessary
notation and preliminary results, as well as formulate the main results.
Section 3 is devoted to the proof of the main result. In Section 4 we discuss
the application of the main results of the paper.

##  2 Preliminaries and Main Results

Throughout this paper $H$ and $K$ are Hibert spaces with norms $\|\cdot\|$ and
$\|\cdot\|_{K}$. Let $(\Omega,\mathcal{F},P)$ be a complete probability space,
and $Q$ be linear bounded covariance operator such that $tr(Q)<\infty.$
Introduce

$W(t):=\sum_{i=1}^{\infty}\sqrt{\lambda_{i}}\beta_{i}(t)e_{i}(x),\ t\geq 0,$

which is a $Q$-Wiener process on $t\geq 0$. Here $\beta_{i}(t)$ are standard,
one dimensional, independent Wiener processes, $\\{e_{k},k\geq 1\\}$ is an
orthonormal system in $K$, and a sequence of nonegative real numbers
$\lambda_{k}$ satisfying

$Qe_{k}=\lambda_{k}e_{k},k=1,2,...$

and

$\sum_{i=1}^{\infty}\lambda_{i}<\infty.$

Also let $\\{F_{t},t\geq 0\\}$ be a normal filtration satisfying

  * •

$W(t)$ is $\mathcal{F}_{t}$-measurable;

  * •

$W(t+h)-W(t)$ is independent of $\mathcal{F}_{t}$ $\forall h\geq 0,t\geq 0$.

As defined in the introduction, let

$\mathcal{L}_{2}^{0}=\mathcal{L}_{2}(Q^{1/2}K,H)$

be the space of all Hilbert Schmidt operators from $Q^{1/2}K$ to $H$ with the
dot product

$\langle{\Phi,\Psi\rangle}_{\mathcal{L}_{2}^{0}}={\rm tr}(\Phi Q\Psi^{*}).$

Let $A$ be an infinitesimal generator of an analytic semigroup $\\{S(t),t\geq
0\\}$ in $H$.

Condition (H1). If $\sigma(-A)\text{ is the spectrum of }(-A)$, we have

${\rm Re}\ \sigma(-A)>\delta>0,\ \text{ and }\ A^{-1}\text{ is compact in }H.$

It follows from [25] that for $0\leq\alpha\leq 1$ one can define the
fractional power $(-A)^{\alpha}$, which is a closed linear operator with
domain $D(-A)^{\alpha}$. Denote $H_{\alpha}$ to be the Banach space
$D(-A)^{\alpha}$ endowed with the norm

$\|u\|_{\alpha}:=\|(-A)^{\alpha}u\|,$

which is equivalent to the graph norm of $(-A)^{\alpha}$. This way $H_{0}=H$.
It follows from [11], Sec. 1.4, that if $A^{-1}$ is compact, then $S(t)$ is
compact for $t>0$. Next, it follows from [25], Th. 3.2, p.48 that under the
assumption (H1) the semigroup $S(t)$ is continuous with respect to uniform
operator topology for $t>0$. Thus, using [25], Th. 3.3, p.48 we may conclude
that the operator $A$ has a compact resolvent. Consequently, using [11], Th.
1.4.8, we have the following result:

######  Lemma 2.1.

Under the condition (H1) the embedding $H_{\alpha}\subset H_{\beta}$ is
compact if $0\leq\beta<\alpha\leq 1$.

######  Lemma 2.2.

([11], Th. 1.4.3) Under the condition (H1), for every $\alpha\geq 0$ there
exists $C_{\alpha}>0$ such that

$\|(-A)^{\alpha}S(t)\|\leq C_{\alpha}t^{-\alpha}e^{-\delta t},t>0.$

In particular,

$\|S(t)\|\leq C_{0}e^{-\delta t},t>0.$ (2)

######  Lemma 2.3.

([5]) Let $p>2$, $T>0$ and let $\Phi$ be an $\mathcal{L}_{2}^{0}$ valued,
predictable process such that

$\mathbb{E}\int_{0}^{T}\|\Phi(t)\|^{p}_{\mathcal{L}_{2}^{0}}dt<\infty.$

Then there is a constant $M_{T}>0$ such that

$\mathbb{E}\sup_{t\in[0,T]}\left\|\int_{0}^{t}S(t-s)\Psi(s)dW(s)\right\|^{p}\leq
M_{T}\mathbb{E}\int_{0}^{T}\|\Psi(s)\|^{p}_{\mathcal{L}_{2}^{0}}ds.$

The following are the main assumptions on the coefficients of the equation
(1):
Condition (H2). The mappings $f:C_{h}\to H$ and
$\sigma:C_{h}\to\mathcal{L}_{2}^{0}$ are continuous and satisfy

  * [i]

There exists a positive constant $K>0$ such that

$\|f(\varphi)\|+\|\sigma(\varphi)\|_{\mathcal{L}_{2}^{0}}\leq
K(1+\|\varphi\|_{C})\text{ for all }\varphi\in C_{h}.$

  * [ii]

There is a function $N:[0,+\infty)\to[0,+\infty)$ such that

    * (a)

the function $N$ is continuous, nondecreasing, concaved, and $N(0)=0$;

    * (b)

For $p>2$ and for all $\varphi_{1},\varphi_{2}\in C_{h}$ we have

$\|f(\varphi_{1})-f(\varphi_{2})\|^{p}+\|\sigma(\varphi_{1})-\sigma(\varphi_{2})\|^{p}_{\mathcal{L}_{2}^{0}}\leq
N\left(\|\varphi_{1}-\varphi_{2}\|^{p}_{C_{h}}\right).$

    * (c)

The only nonnegative function $z(t)$ satisfying

$z(t)\leq D\int_{0}^{t}N(z(s))ds\text{ for all }t\in[0,T]$

is $z\equiv 0$.

######  Remark 2.4.

The condition (c) holds for any $N$ satisfying (a) and

$\int_{0}^{1}\frac{ds}{N(s)}=+\infty.$

Condition (H3). There exist positive constants $0<\alpha\leq 1$ and
$0<M_{g}<1$ such that for all $\varphi_{1},\varphi_{2}\in C_{h}$ the function
$g:C_{h}\to H_{\alpha}$ satisfies

$\|g(\varphi_{1})-g(\varphi_{2})\|_{H_{\alpha}}\leq
M_{g}\|\varphi_{1}-\varphi_{2}\|_{C_{h}}.$

Condition (H4). The initial condition $\varphi:[-h,0]\times\Omega\to H$ is an
$\mathcal{F}_{0}$ \- measurable random variable, which is independent of $W$
and has continuous trajectories.

######  Definition 2.5.

[18, 13, 8] A continuous $\mathcal{F}_{t}$-adapted stochastic process
$u:[-h,T]\times\Omega\to H$ is a mild solution for (1) for $t\in[0,T]$ if it
satisfies the integral equation

$u(t)=S(t)(\varphi(0)+g(\varphi))-g(u_{t})-\int_{0}^{t}AS(t-s)g(u_{s})ds+\int_{0}^{t}S(t-s)f(u_{s})ds+\int_{0}^{t}S(t-s)\sigma(u_{s})dW(s)$
(3)

and $u(t)=\varphi(t)$ a.s. for $t\in[-h,0]$.

The main results of this paper are the Theorem 2.6, Theorem 2.8 and Theorem
2.13 below.

######  Theorem 2.6.

(Existence and uniqueness) Suppose the conditions (H1)-(H4) hold. Then for all
$T>0$ the equation (1) has a unique mild solution $u$ on $[0,T]$. Furthermore,
this solution satisfies

$\|u_{t}-u_{0}\|_{C_{h}}=\sup_{\theta\in[-h,0]}\|u(t+\theta)-\varphi(\theta)\|\to^{P}0,t\to
0.$

The main idea of the proof of Theorem 2.6 is to consider an auxiliary
semilinear equation, which, in turn, will enable us to constuct a convergent
approximating sequence. In contrast with non-neutral type equations (with
$g\equiv 0$), the primary difficulty of the equation (1) comes from the term
$AS(t-s)g(u_{s})$. Generally speaking, if $g(u_{s})\in H$, this term has a
non-integrable singularity at $t=s$. We overcome this difficulty by
introducing the fractional powers of the operator $(-A)$. In particular, we
show that if $g$ is regular enough for $g(u_{s})\in D(-A)^{\alpha}$ for
$\alpha>0$, then this singularity becomes integrable.

######  Remark 2.7.

In order to prove Theorem 2.6 we will use the methods outlined in [27], which
established the corresponding result for the neutral type stochastic equations
in a simpler case. Besides, Theorem 2.6 has a different condition on the
Lipschitz constant for $g$. Namely, [27] requires

$4^{p-1}\|(-A)^{-\alpha}\|^{p}M_{g}<1,$

while we require a much more easy to check condition $M_{g}<1$. Similar
conditions on $(-A)^{-\alpha}$ appear in [13, 1] among others. In addition, we
relax the regularity condition for the map $g(\varphi)$. In particular, we
require $g\in H_{\alpha}$ for $\alpha\in(0,1)$, while in [27] the authors
require $\alpha\in(\frac{1}{p},1)$ for $p>2$.

######  Theorem 2.8.

(Continuous dependence on the initial data) Suppose the conditions of Theorem
2.6 hold. Let $u(t,\varphi)$ and $u(t,\Psi)$ be two solutions of (1) with the
initial conditions $\varphi$ and $\Psi$ satisfying the condition (H4). Then

$\mathbb{E}\sup_{t\in[0,T]}\left[\|u(t,\varphi)-u(t,\Psi)\|^{p}+\|u_{t}(\varphi)-u_{t}(\Psi)\|_{C_{h}}^{p}\right]\to
0,$

as $\mathbb{E}\|\varphi-\Psi\|_{C_{h}}^{p}\to 0.$

Denote $B_{b}(C_{h})$ to be the Banach space of bounded real Borel measurable
functionals, on $C_{h}$, and $C_{b}(C_{h})$ is the space of bounded continuous
functionals on $C_{h}$. Since Theorem 2.6 guarantees the existence and
uniqueness of the solution at any $t\geq 0$, replacing the initial interval
$[-h,0]$ with $[-h+s,s]$ for any $s\geq 0$ we can guarantee the existence and
uniqueness of the solution on any interval $[s,t]$ with the initial
$\mathcal{F}_{s}$ measurable function $\varphi$ satisfying the conditions on
$[s-h,s]$. Such solutions will be denoted with $u(t,s,\varphi)$. Similarly,
denote $u_{t}(s,\theta)=u(t+\theta,s,\varphi),\theta\in[-h,0]$ to be the shift
of the solution $u(t,\varphi)$, such that
$u_{s}(s,\varphi)=u(s+\theta,s,\varphi)=\varphi(\theta)$.

Following [31], define the family of shift operators

$U_{s}^{t}\varphi:=u(t+\theta,s,\varphi)=u_{t}(s,\varphi).$ (4)

Denote $\mathcal{F}_{s}^{t}(dW)$ to be the minimal $\sigma$ \- algebra
containing $W(\tau)-W(s),\tau\in[s,t]$, and $G^{t}$ to be the minimal $\sigma$
\- algebra containing $W(\tau)-W(t)$ for $\tau\geq t$. For any nonrandom
$\varphi\in C_{h}$ with $s\geq 0$ and $t\geq s$ the solution
$U_{s}^{t}\varphi=u_{t}(s,\varphi)$ is a $\mathcal{F}_{s}^{t}(dW)$\-
measurable random function taking values in $C_{h}$. Note that in this case,
since $\varphi$ is nonrandom, $u_{t}(s,\varphi)$ is independent of the
$\sigma$ \- algebra $G^{t}$. The next proposition was shown, e.g. in [21]:

######  Proposition 2.9.

The operator (4) satisfies

$U_{\tau}^{t}U_{s}^{\tau}\varphi=U_{s}^{t}\varphi$

for all $t\geq\tau\geq s\geq 0$, and $\varphi\in C_{h}$.

Let $D$ be a $\sigma$-algebra of Borel subsets of $C_{h}$. For every $A\in D$
we define

$\mu_{t}(A)=P\\{u_{t}(s,\varphi)\in A\\}=P\\{U_{s}^{t}\varphi\in
A\\}=P(s,\varphi,t,A).$ (5)

This way $u_{t}(s,\varphi)$ naturally defines a measure on $D$. The formula
(5) defines a transition probability function, corresponding to the random
process $u_{t}(s,\varphi),t\geq s\geq 0$. In a similar way to the finite
dimensional case [31], it is possible to show that this function satisfies the
properties of the transition probability. This way we have

######  Theorem 2.10.

(Markov property) Under the assumptions of Theorem 2.6, the process
$u_{t}(s,\varphi)$ is the Markov process on $C_{h}$ with the transition
function, defined by (5).

######  Proposition 2.11.

For any $t\geq s\geq 0$ and any $A\in D$ we have

$P(s,\varphi,t,A)=P(0,\varphi,t-s,A).$

###### Proof.

Let $\tilde{u}(t)=u(s+t,s,\varphi)$. Then
$\tilde{u}_{0}=u(s+\theta,s,\varphi),\tilde{u}_{t}=u(s+t+\theta,s,\varphi)=u_{s+t}(s,\varphi)$,
for $\theta\in[-h,0].$ On the other hand,

$\displaystyle\tilde{u}(t)=u(s+t,s,\varphi)=S(t)(\varphi(0)+g(\varphi))-g(u_{t+s})-\int_{s}^{s+t}AS(t+s-\tau)g(u_{\tau})d\tau$
$\displaystyle+\int_{s}^{s+t}S(t+s-\tau)f(u_{\tau})d\tau+\int_{s}^{s+t}S(t+s-\tau)\sigma(u_{\tau})dW(\tau)=S(t)(\varphi(0)+g(\varphi))-g(u_{t+s})$
$\displaystyle-\int_{0}^{t}AS(t-\tau)g(u_{\tau+s})d\tau+\int_{0}^{t}S(t-\tau)f(u_{\tau+s})d\tau+\int_{0}^{t}S(t-\tau)\sigma(u_{\tau+s})d\tilde{W}(\tau),$

where $\tilde{W}(\tau):=W(s+\tau)-W(s)$ is a $Q$-Wiener process as well. This
way

$\displaystyle\tilde{u}(t)=S(t)(\varphi(0)+g(\varphi))-g(\tilde{u}_{t})-\int_{0}^{t}AS(t-\tau)g(\tilde{u}_{\tau})d\tau$
(6)
$\displaystyle+\int_{0}^{t}S(t-\tau)f(\tilde{u}_{\tau})d\tau+\int_{0}^{t}S(t-\tau)\sigma(\tilde{u}_{\tau})d\tilde{W}(\tau).$

However, the same equation is satisfied with $u(t,0,\varphi)$, with
$u(0,0,\varphi)=\varphi(0)$ and $u_{0}=\varphi(\theta)$. The only difference
is that $u(t,0,\varphi)$ solves (6) with a different $Q$-Wiener process
$\tilde{W}$. However, since the distribution function for $W$ is the same as
for $\tilde{W}$, the uniqueness of solution implies that the distributions of
$u(s+t,s,\varphi)$ and $u(t,0,\varphi)$ coincide. In other words, the
distribution of $u(s+t,s,\varphi)$ is independent of $s$. Then the
distribution of
$u_{t}(s,\varphi)=u(t+\theta,s,\varphi)=u(t-s+\theta+s,s,\varphi)$ matches the
distribution of $u(t-s+\theta,0,\varphi)=u_{t-s}(0,\varphi)$. Hence

$\displaystyle P(s,\varphi,t,A)=P\\{u_{t}(s,\varphi)\in A)\\}$
$\displaystyle=$ $\displaystyle P\\{u(t+\theta,s,\varphi)\in A\\}$
$\displaystyle=$ $\displaystyle P\\{u(t+\theta-s,0,\varphi)\in
A\\}=P\\{u_{t-s}(0,\varphi)\in A\\},$

which yields the desired result. ∎

For $g\in B_{b}(C_{h})$, for all $\varphi\in C_{h}$ and $t\geq s\geq 0$ define

$P_{s,t}(\varphi):=\mathbb{E}\,g(u_{t}(s,\varphi)).$

From Proposition 2.11 we have $P_{s,t}(\varphi)=P_{0,t-s}(\varphi)$ and denote
$P_{t}\varphi=P_{0,t}(\varphi)$. The following proposition follows from
Theorem 2.8.

######  Proposition 2.12.

Under the assumptions (H1)-(H4) the transition semigroup $P_{t},t\geq 0$ is
stochastically continuous and and satisfies the Feller property

$P_{t}:C_{b}(C_{h})\to C_{b}(C_{h})\text{ and }\lim_{t\to
0}P_{t}\varphi=\varphi.$

We are now in position to state our main result on the long time behavior of
solutions of (1).

######  Theorem 2.13.

Assume the assumptions of Theorem 2.6 hold. Suppose the equation (1) has a
solution $u(t)$ which is bounded in probability for $t\geq 0$ in $C_{h}$,
namely

$\sup_{t\geq 0}P\\{\|u_{t}\|_{C_{h}}>R\\}\to 0,R\to\infty.$ (7)

Then there exists an invariant measure $\mu$ in $C_{h}$, i.e.

$\int_{C_{h}}P_{t}\varphi(x)d\mu(x)=\int_{C_{h}}\varphi(x)d\mu(x),\text{ for
any }t\geq 0,\text{ and }\varphi\in C_{b}(C_{h}).$

######  Remark 2.14.

The condition (7) is a standard condition for the existence of invariant
measures, see e.g. [4]. In Section 4 we will show how to verify this condition
in terms of the coefficients of a neutral type stochastic functional-
differential equation with elliptic operator.

##  3 Proofs of the main results

###  3.1 Proof of existence and uniqueness

In this subsection we will prove Theorem 2.6. We are going to use the approach
from [27]. Assume first that
$\mathbb{E}\sup_{\theta\in[-h,0]}\|\varphi(\theta)\|^{p}<\infty,p>2$. Let
$\mathcal{L}_{\mathcal{F}}^{p}([0,T],H)$ be the Banach space of $H$-valued,
$\mathcal{F}_{t}$-adapted random processes with the norm

$\|\zeta\|^{p}_{H,T}:=\mathbb{E}\int_{0}^{T}\|\zeta(t)\|^{p}dt.$

Similarly, let $\mathcal{L}_{\mathcal{F}}^{p}([0,T],\mathcal{L}_{2}^{0})$, be
the Banach space of $\mathcal{L}_{2}^{0}$ valued $\mathcal{F}_{t}$-adapted
random processes $\Phi(t)$ with the norm

$\|\Phi\|^{p}_{L_{2}^{0},T}:=\mathbb{E}\int_{0}^{T}\|\Phi(t)\|_{L_{2}^{0}}^{p}dt.$

Next, denote $\mathcal{B}_{p,T}$ to be the Banach space of all $H$-valued,
measurable, continuous $\mathcal{F}_{t}$-adapted (for $t\geq 0$) processes
$u(t,\omega):[-h,T]\times\Omega\to H$, equipped with the norm

$\|u\|_{\mathcal{B}_{p,T}}^{p}:=\mathbb{E}\sup_{t\in[-h,T]}\|u(t)\|^{p}.$

Finally, introduce the closed set

$\mathcal{B}_{p,T}(\varphi):=\\{u\in\mathcal{B}_{p,T},u(t)=\varphi(t)\text{
for }t\in[-h,0]\\}.$

We will use the following Lemma:

######  Lemma 3.1.

For $T>0$, let $u(t),-h\leq t\leq T$ be a stochastic process with continuous
paths. Then, if $p\geq 1$ and $u_{0}=\varphi$, we have

$\mathbb{E}\sup_{t\in[0,T]}\|u_{t}\|_{C_{h}}^{p}\leq\mathbb{E}\|\varphi\|_{C_{h}}^{p}+\mathbb{E}\sup_{t\in[0,T]}\|u(t)\|^{p}.$

The proof of Lemma 3.1 follows directly from the definition of $u_{t}$. We
next split the proof into the following steps:
Step 1. We start with an auxiliary equation. For any
$(a,b)\in\mathcal{L}_{\mathcal{F}}^{p}([0,T],H)\times\mathcal{L}_{\mathcal{F}}^{p}([0,T],L_{2}^{0})$
consider the following equation

$\begin{cases}d(u(t)+g(u_{t}))=[Au(t)+a(t)]dt+b(t)dW(t),t\geq 0;\\\
u(t)=\varphi(t),t\in[-h,0].\end{cases}$ (8)

######  Lemma 3.2.

Under the condition (H3) the equation (8) has a unique mild solution in
$\mathcal{B}_{p,T}(\varphi)$.

###### Proof.

Define the following operator
$\Psi:\mathcal{B}_{p,T}(\varphi)\to\mathcal{B}_{p,T}(\varphi)$ as

$[\Psi
u](t):=S(t)(\varphi(0)+g(\varphi))-g(u_{t})-\int_{0}^{t}AS(t-s)g(u_{s})ds+\int_{0}^{t}S(t-s)a(s)ds+\int_{0}^{t}S(t-s)b(s)dW(s),$
(9)

defined for $t\in[0,T]$, with $u(t)=\varphi(t)$ for $t\in[-h,0]$. Let us show
that $\Psi:\mathcal{B}_{p,T}(\varphi)\to\mathcal{B}_{p,T}(\varphi)$ is a
contraction. To this end, let us first show that $\Psi$ takes values in
$\mathcal{B}_{p,T}(\varphi)$. Using Lemma 3.1 we have

$\mathbb{E}\sup_{t\in[0,T]}\|g(u_{t})\|^{p}\leq\mathbb{E}\sup_{t\in[0,T]}\|g(u_{t})\|^{p}_{\alpha}\leq
M_{g}^{p}(\|g(0)\|_{C_{h}}^{p}+\mathbb{E}\sup_{t\in[0,T]}\|u(t)\|^{p}).$

Next,

$\displaystyle\mathbb{E}\sup_{t\in[0,T]}\left\|\int_{0}^{t}AS(t-s)g(u_{s})ds\right\|^{p}=\mathbb{E}\sup_{t\in[0,T]}\left\|\int_{0}^{t}(-A)^{1-\alpha}S(t-s)(-A)^{\alpha}g(u_{s})ds\right\|^{p}$
$\displaystyle\leq\mathbb{E}\sup_{t\in[0,T]}\left(\int_{0}^{t}\left\|(-A)^{1-\alpha}S(t-s)\right\|\left\|(-A)^{\alpha}g(u_{s})\right\|ds\right)^{p}$
$\displaystyle\leq\sup_{t\in[0,T]}\left(\int_{0}^{t}C_{1-\alpha}(t-s)^{-(1-\alpha)}ds\right)^{p}\mathbb{E}\sup_{t\in[0,T]}\|g(u_{t})\|^{p}_{\alpha}.$

The remaining two terms can be estimated in a standard way (see, e.g. [27]):

$\mathbb{E}\sup_{t\in[0,T]}\left\|\int_{0}^{t}S(t-s)a(s)ds\right\|^{p}\leq
C_{0}^{p}\|a\|^{p}_{H,T}$

and

$\mathbb{E}\sup_{t\in[0,T]}\left\|\int_{0}^{t}S(t-s)b(s)dW(s)\right\|^{p}\leq
C_{0}^{p}\|b\|^{p}_{L_{2}^{0},T},$

where $C_{0}$ is given by (2).

Let us now show that $\Psi$ is contraction. Fix $T_{1}\in[0,T]$ For any
$\rho\in(0,1)$, and $a,b\geq 0$ we will use the inequality

$(a+b)^{p}\leq\frac{a^{p}}{\rho^{p-1}}+\frac{b^{p}}{(1-\rho)^{p-1}}.$

For any $u,v\in B_{p,T_{1}}$, we have

$\displaystyle\left\|\Psi(u)-\Psi(v)\right\|^{p}_{B_{p,T_{1}}}\leq\frac{1}{M_{g}^{p-1}}\|g(u_{t})-g(v_{t})\|^{p}_{B_{p,T_{1}}}$
$\displaystyle+\frac{1}{(1-M_{g})^{p-1}}\left\|\int_{0}^{t}AS(t-s)(g(u_{s})-g(v_{s}))ds\right\|^{p}_{B_{p,T_{1}}}\leq\frac{1}{M_{g}^{p-1}}\mathbb{E}\sup_{t\in[0,T_{1}]}\|g(u_{t})-g(v_{t})\|^{p}$
$\displaystyle+\frac{1}{(1-M_{g})^{p-1}}\mathbb{E}\sup_{t\in[0,T_{1}]}\left(\int_{0}^{t}\|(-A)^{1-\alpha}S(t-s)\|\|(-A)^{\alpha}(g(u_{s})-g(v_{s})\|ds\right)^{p}$
$\displaystyle\leq
M_{g}\mathbb{E}\sup_{t\in[0,T_{1}]}\|u_{t}-v_{t}\|^{p}_{C_{h}}+\frac{M_{g}^{p}C_{1-\alpha}^{p}T_{1}^{\alpha
p}}{(1-M_{g})^{p-1}\alpha^{p}}\mathbb{E}\sup_{t\in[0,T_{1}]}\|u_{t}-v_{t}\|^{p}_{C_{h}}=\gamma(T_{1})\mathbb{E}\sup_{t\in[0,T_{1}]}\|u_{t}-v_{t}\|^{p}_{C_{h}}.$

Here we used Lemma 3.1 and Lemma 2.2. Using the assumption (H3) we can choose
$T_{1}$ such that

$\gamma(T_{1}):=M_{g}+\frac{M_{g}^{p}C_{1-\alpha}^{p}T_{1}^{\alpha
p}}{(1-M_{g})^{p-1}\alpha^{p}}<1,$ (10)

hence the operator $\Psi$ is a contraction in
$\mathcal{B}_{p,T_{1}}(\varphi)$. Therefore $\Psi$ has a unique fixed point,
which is a unique fixed point, which is a mild solution of (8) on $[0,T_{1}]$.

Let us show that the solution has a continuous trajectory. It is
straightforward to verify the continuity of the first, third and fourth terms
in the equation (9). The continuity of the fifth term follows from the
factorization method [5], Prop. 7.3.

Let us now establish the continuity of the second term. Indeed, since $u$ is
continuous,

$\|g(u_{t+s})-g(u_{t})\|\leq
M_{g}\|u_{t+s}-u_{t}\|_{C_{h}}=\sup_{\theta\in[-h,0]}\|u(t+s+\theta)-u(t+\theta)\|\to
0,s\to 0\ a.s.$

We may now consider the new initial value problem with the initial condition
on $[T_{1}-h,T_{1}]$ and so on. This way the solution may be extended to the
entire interval $[0,T]$ in finitely many steps. Hence the equation (8) has a
unique mild solution, defined on $[0,T]$, which completes the proof of the
Lemma. ∎

Step 2. For given $\varphi\in C_{h}$ we define the operator

$\Phi_{\varphi}:L^{p}_{\mathcal{F}}([0,T],H)\times
L^{p}_{\mathcal{F}}([0,T],L_{2}^{0})\to\mathcal{B}_{p,T}$

as

$(\Phi_{\varphi}(a,b))(t):=\begin{cases}\varphi(t),t\in[-h,0];\\\
u(t),t\in[0,T],\end{cases}$

where $u$ is the unique solution of (8).

######  Lemma 3.3.

Under the assumption (H3), for every $p>2$, there are positive constants $B$,
$D$ and $L$, depending only on $T$, $M_{g}$, $p$ and $\alpha$, such that for
any $t\in[0,T]$, for any $\varphi$ and $\psi$ in $C_{h}$, and for all pairs

$(a,b),(a_{1},b_{1})\text{ and }(a_{2},b_{2})\in
L^{p}_{\mathcal{F}}([0,T],H)\times L^{p}_{\mathcal{F}}([0,T],L_{2}^{0}),$

we have

$\displaystyle\|\Phi_{\varphi}(a_{1},b_{1})-\Phi_{\psi}(a_{2},b_{2})\|^{p}_{\mathcal{B}_{p,t}}\leq
L\|\varphi-\psi\|^{p}_{C_{h}}$ (11)
$\displaystyle+B\int_{0}^{t}\left[\mathbb{E}\|a_{1}(s)-a_{2}(s)\|^{p}+\mathbb{E}\|b_{1}(s)-b_{2}(s)\|^{p}_{L_{2}^{0}}\right]ds,\
\ \text{ and }$ $\displaystyle\|\Phi(a,b)\|^{p}_{\mathcal{B}_{p,t}}\leq
D+B\int_{0}^{t}\left[\mathbb{E}\|a(s)\|^{p}+\mathbb{E}\|b(s)\|^{p}_{L_{2}^{0}}\right]ds.$
(12)

###### Proof.

For any $(a_{1},b_{1})$ and $(a_{2},b_{2})$ we have

$\Phi_{\varphi}(a_{1},b_{1}):=u=S(t)(\varphi(0)+g(\varphi))-g(u_{t})-\int_{0}^{t}AS(t-s)g(u_{s})ds\\\
+\int_{0}^{t}S(t-s)a_{1}(s)ds+\int_{0}^{t}S(t-s)b_{1}(s)dW(s),$

and

$\Phi_{\psi}(a_{2},b_{2}):=v=S(t)(\psi(0)+g(\psi))-g(v_{t})-\int_{0}^{t}AS(t-s)g(v_{s})ds\\\
+\int_{0}^{t}S(t-s)a_{2}(s)ds+\int_{0}^{t}S(t-s)b_{2}(s)dW(s).$

Denote

$\beta(s):=\mathbb{E}\|a_{1}(s)-a_{2}(s)\|^{p}+\mathbb{E}\|b_{1}(s)-b_{2}(s)\|^{p}_{L_{2}^{0}}.$

Now choose $T_{1}>0$, which satisfies both (10) and

$M_{g}+\left(\frac{5}{1-M_{g}}\right)^{p-1}\left(\frac{C_{1-\alpha}T_{1}^{\alpha}M_{g}}{\alpha}\right)^{p}<1.$
(13)

Note that such $T_{1}$ always exists provided $M_{g}<1$. Next, we partition
the interval $[0,T]$ with the points $kT_{1}$, and consider the interval
$t\in[0,T_{1}]$ first. For such $t$ we have

$\displaystyle\|u-v\|_{\mathcal{B}_{p,t}}^{p}=\mathbb{E}\sup_{s\in[0,t]}\|u(s)-v(s)\|^{p}\leq\mathbb{E}\sup_{s\in[0,t]}\Big{[}\|S(s)(g(\varphi)-g(\psi))\|+\|S(s)(\varphi(0)-\psi(0))\|$
$\displaystyle+\|g(u_{s})-g(v_{s})\|+\left\|\int_{0}^{s}AS(s-\tau)(g(u_{\tau})-g(v_{\tau}))d\tau\right\|+\left\|\int_{0}^{s}S(s-\tau)(a_{1}(\tau)-a_{2}(\tau))d\tau\right\|$
$\displaystyle+\left\|\int_{0}^{s}S(s-\tau)(b_{1}(\tau)-b_{2}(\tau))dW(\tau)\right\|\Big{]}^{p}\leq\frac{1}{M_{g}^{p-1}}\mathbb{E}\sup_{s\in[0,t]}\|g(u_{s})-g(v_{s})\|^{p}$
(14)
$\displaystyle+\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}(1+M_{g}^{p})\|\varphi-\psi\|_{C_{h}}^{p}$
$\displaystyle+\left(\frac{5}{1-M_{g}}\right)^{p-1}\mathbb{E}\sup_{s\in[0,t]}\left(\int_{0}^{s}\|(-A)^{1-\alpha}S(s-\tau)\|\|(g(u_{\tau})-g(v_{\tau})\|_{\alpha}d\tau\right)^{p}$
$\displaystyle+C_{0}^{p}\left(\frac{5}{1-M_{g}}\right)^{p-1}\int_{0}^{t}\beta(s)ds,$

where $C_{0}$ is given in (2). Let us estimate each term in (3.1) separately.

$\displaystyle\mathbb{E}\sup_{s\in[0,t]}\|g(u_{s})-g(v_{s})\|^{p}\leq\mathbb{E}\sup_{s\in[0,t]}\|g(u_{s})-g(v_{s})\|^{p}_{\alpha}\leq
M_{g}^{p}\mathbb{E}\sup_{s\in[0,t]}\|u_{s}-v_{s}\|^{p}_{C_{h}}$
$\displaystyle=M_{g}^{p}\mathbb{E}\sup_{s\in[0,t]}\sup_{\theta\in[-h,0]}\|u(s+\theta)-v(s+\theta))\|^{p}=M_{g}^{p}\mathbb{E}\sup_{s\in[-h,t]}\|u(s)-v(s)\|^{p}$
$\displaystyle\leq
M_{g}^{p}\|\varphi-\psi\|_{C_{h}}^{p}+M_{g}^{p}\mathbb{E}\sup_{s\in[0,t]}\|u(s)-v(s)\|^{p}.$

Next,

$\mathbb{E}\sup_{s\in[0,t]}\left(\int_{0}^{s}\|(-A)^{1-\alpha}S(s-\tau)\|\|(g(u_{\tau})-g(v_{\tau})\|_{\alpha}d\tau\right)^{p}\leq\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}\mathbb{E}\sup_{s\in[-h,t]}\|u(s)-v(s)\|^{p}\\\
\leq\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}\left(\mathbb{E}\sup_{s\in[0,t]}\|u(s)-v(s)\|^{p}+\|\varphi-\psi\|_{C_{h}}^{p}\right).$
(15)

Using (3.1) - (15), we have

$\|u-v\|^{p}_{\mathcal{B}_{p,t}}\leq
M_{g}\mathbb{E}\sup_{s\in[0,t]}\|u(s)-v(s)\|^{p}+\left(\frac{5}{1-M_{g}}\right)^{p-1}\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}\mathbb{E}\sup_{s\in[0,t]}\|u(s)-v(s)\|^{p}\\\
+\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}\int_{0}^{t}\beta(s)ds+\tilde{L}_{1}\|\varphi-\psi\|^{p}_{C_{h}},$

or

$\|u-v\|^{p}_{\mathcal{B}_{p,t}}\leq
C_{T_{1}}\int_{0}^{t}\beta(s)ds+L_{1}\|\varphi-\psi\|^{p}_{C_{h}},$ (16)

where

$C_{T_{1}}=\frac{\left(\frac{5}{1-M_{g}}\right)^{p-1}}{1-M_{g}-\left(\frac{5}{1-M_{g}}\right)^{p-1}\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}}C_{0}^{p}$

and

$L_{1}=\left(\frac{5}{1-M_{g}}\right)^{p-1}\frac{C_{0}^{p}+M_{g}^{p}C_{0}^{p}+M_{g}^{p}+\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}}{1-M_{g}-\left(\frac{5}{1-M_{g}}\right)^{p-1}\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}}.$

Let us now consider the second interval $t\in[T_{1},2T_{1}]$. We have

$u(t)=S(t-T_{1})(u(T_{1})+g(u_{T_{1}}))-g(u_{t})-\int_{T_{1}}^{t}AS(t-s)g(u_{s})ds\\\
+\int_{T_{1}}^{t}S(t-s)a_{1}(s)ds+\int_{T_{1}}^{t}S(t-s)b_{1}(s)dW(s),$ (17)

and

$v(t)=S(t-T_{1})(v(T_{1})+g(v_{T_{1}}))-g(v_{t})-\int_{T_{1}}^{t}AS(t-s)g(v_{t})ds\\\
+\int_{T_{1}}^{t}S(t-s)a_{2}(s)ds+\int_{T_{1}}^{t}S(t-s)b_{2}(s)dW(s).$

Then

$\displaystyle\|u-v\|_{\mathcal{B}_{p,t}}^{p}=$
$\displaystyle\mathbb{E}\sup_{s\in[0,t]}\|u(s)-v(s)\|^{p}$
$\displaystyle=\max\\{\mathbb{E}\sup_{s\in[0,T_{1}]}\|u(s)-v(s)\|^{p},\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}\\}$
$\displaystyle\leq\mathbb{E}\sup_{s\in[0,T_{1}]}\|u(s)-v(s)\|^{p}+\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}.$

Let us now estimate the second term in (3.1). From (17) we have

$\displaystyle\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}\leq\mathbb{E}\sup_{s\in[T_{1},t]}\Big{(}\|g(u_{s})-g(v_{s})\|+\|S(s-T_{1})\|(\|u(T_{1})-v(T_{1})\|$
$\displaystyle+\|g(u_{T_{1}})-g(v_{T_{1}})\|)+\int_{T_{1}}^{s}C_{1-\alpha}(s-\tau)^{-(1-\alpha)}M_{g}\|u_{\tau}-v_{\tau}\|_{C_{h}}d\tau+$
$\displaystyle\left\|\int_{T_{1}}^{s}S(s-\tau)(a_{1}(\tau)-a_{2}(\tau))d\tau\right\|+\left\|\int_{T_{1}}^{s}S(s-\tau)(b_{1}(\tau)-b_{2}(\tau))dW(\tau)\right\|\Bigg{)}^{p}$
(19)
$\displaystyle\leq\frac{1}{M_{g}^{p-1}}\mathbb{E}\sup_{s\in[T_{1},t]}\|g(u_{s})-g(v_{s})\|^{p}+\frac{5^{p-1}}{(1-M_{g})^{p-1}}\Big{(}C_{0}^{p}\Big{(}\mathbb{E}\|u(T_{1})-v(T_{1})\|^{p}$
$\displaystyle\mathbb{E}\|g(u_{T_{1}})-g(v_{T_{1}})\|^{p}\Big{)}\Big{)}+\mathbb{E}\sup_{s\in[T_{1},t]}\left(\int_{T_{1}}^{s}C_{1-\alpha}(s-\tau)^{-(1-\alpha)}M_{g}\|u_{\tau}-v_{\tau}\|_{C_{h}}d\tau\right)^{p}+C_{0}^{p}\int_{T_{1}}^{t}\beta(\tau)d\tau.$

Setting $t=T_{1}$ in (16), we have

$\mathbb{E}\|u(T_{1})-v(T_{1})\|^{p}\leq
C_{T_{1}}\int_{0}^{T_{1}}\beta(s)ds+L_{1}\|\varphi-\psi\|_{C_{h}}^{p}$

and

$\mathbb{E}\|g(u_{T_{1}})-g(v_{T_{1}})\|^{p}\leq
M_{g}^{p}\mathbb{E}\sup_{\theta\in[-h,0]}\|u(T_{1}+\theta)-v(T_{1}+\theta)\|^{p}\\\
\leq M_{g}^{p}\mathbb{E}\sup_{s\in[-h,T_{1}]}\|u(s)-v(s)\|^{p}\leq
M_{g}^{p}C_{T_{1}}\int_{0}^{T_{1}}\beta(s)ds+L_{21}\|\varphi-\psi\|_{C_{h}}^{p}.$

Using the result of Lemma 3.1, we proceed with the following estimate:

$\displaystyle\mathbb{E}\sup_{s\in[T_{1},t]}\|g(u_{s})-g(v_{s})\|^{p}\leq
M_{g}^{p}\mathbb{E}\sup_{s\in[T_{1},t]}\sup_{\theta\in[-h,0]}\|u(s+\theta)-v(s+\theta)\|^{p}$
$\displaystyle\leq
M_{g}^{p}\mathbb{E}\sup_{s\in[T_{1}-h,t]}\|u(s)-v(s)\|^{p}\leq
M_{g}^{p}\max\\{\mathbb{E}\sup_{s\in[0,T_{1}]}\|u(s)-v(s)\|^{p},\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}\\}$
$\displaystyle\leq
M_{g}^{p}C_{T_{1}}\int_{0}^{T_{1}}\beta(s)ds+M_{g}^{p}\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}+M_{g}^{p}L_{1}\|\varphi-\psi\|_{C_{h}}^{p}.$

We next estimate the integral term in (3.1). Using Lemma 3.1 once again, we
have

$\mathbb{E}\sup_{s\in[T_{1},t]}\left(\int_{T_{1}}^{s}\frac{C_{1-\alpha}}{(s-\tau)^{1-\alpha}}M_{g}\|u_{\tau}-v_{\tau}\|_{C_{h}}d\tau\right)^{p}$
$\leq\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}\
\mathbb{E}\sup_{s\in[T_{1},t]}(\sup_{\tau\in[T_{1},s]}(\sup_{\theta\in[-h,0]}\|u(\tau+\theta)-v(\tau+\theta)\|^{p}))$
$=\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}\mathbb{E}\sup_{s\in[T_{1}-h,t]}\|u(s)-v(s)\|^{p}\leq\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}C_{T_{1}}\int_{0}^{T_{1}}\beta(s)ds$
(20) $+\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}\
\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}+L_{22}\|\varphi-\psi\|^{p}_{C_{h}}.$

It follows from (3.1)-(20) that

$\displaystyle\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}\leq
M_{g}C_{T_{1}}\int_{0}^{T_{1}}\beta(s)ds+M_{g}\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}$
$\displaystyle\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}C_{T_{1}}\int_{0}^{T_{1}}\beta(s)ds+\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}C_{T_{1}}M_{g}^{p}\int_{0}^{T_{1}}\beta(s)ds$
$\displaystyle+\left(\frac{5}{1-M_{g}}\right)^{p-1}\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}$
$\displaystyle+\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}\int_{T_{1}}^{t}\beta(s)ds+L_{23}\|\varphi-\psi\|_{C_{h}}^{p},$

or,

$\displaystyle\mathbb{E}\sup_{s\in[T_{1},t]}\|u(s)-v(s)\|^{p}\left(1-M_{g}-\left(\frac{5}{1-M_{g}}\right)^{p-1}\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}\right)$
$\displaystyle\leq\left(M_{g}C_{T_{1}}+\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}C_{T_{1}}+\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}C_{T_{1}}M_{g}^{p}\right)\int_{0}^{T_{1}}\beta(s)ds$
$\displaystyle+\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}\int_{T_{1}}^{t}\beta(s)ds+L_{23}\|\varphi-\psi\|_{C_{h}}^{p}.$

In view of (3.1), we conclude that

$\|u-v\|_{\mathcal{B}_{p,t}}^{p}\leq
C_{2}(T_{1})\int_{0}^{t}\beta(s)ds+L_{2}(T_{1})\|\varphi-\psi\|_{C_{h}}^{p},$

where

$C_{2}(T_{1})=\frac{\tilde{C_{2}}(T_{1})}{1-M_{g}-\left(\frac{5}{1-M_{g}}\right)^{p-1}\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}}$

with

$\tilde{C}_{2}(T_{1})=C_{T_{1}}M_{g}+\left(\frac{5}{1-M_{g}}\right)^{p-1}C_{0}^{p}[C_{T_{1}}+C_{T_{1}}M_{g}^{p}+1].$

Repeating the arguments on $[kT_{1},(k+1)T_{1}]$, we have

$\|u-v\|_{\mathcal{B}_{p,t}}^{p}\leq
C_{k}(T_{1})\int_{0}^{t}\beta(s)ds+L_{k}(T_{1})\|\varphi-\psi\|_{C_{h}}^{p},$

where

$C_{k}(T_{1})=\frac{\tilde{C_{k}}(T_{1})}{1-M_{g}-\left(\frac{5}{1-M_{g}}\right)^{p-1}\left(\frac{M_{g}C_{1-\alpha}T_{1}^{\alpha}}{\alpha}\right)^{p}}$

and $\tilde{C}_{k}(T_{1})$ depends only on $C_{i}(T_{1}),i=1,...,k-1$, and,
therefore, is well defined for $T_{1}$ satisfying the conditions (10) and
(13). Similarly, $L_{k}(T_{1})$ depends only on $L_{i}(T_{1}),i=1,...,k-1$,
and, therefore, is also well defined for $T_{1}$ satisfying the conditions
(10) and (13). Since the interval $[0,T]$ consists of $N$ intevals of length
$T_{1}$ (with $T_{1}$ dependent only on $M_{g}$), (11) follows with
$B=\max_{0\leq k\leq N-1}C_{k}$. Finally, the inequality (12) follows from
(11) by setting $a_{2}=b_{2}=0$. ∎

Construction of the approximation sequence. The existence and uniqueness of
mild solution of (1) will be established using Picard type iteration scheme.
For fixed $T>0$, let $u^{(0)}$ be a solution of (8) with $a\equiv b\equiv 0.$
We may now define

$u^{(n)}(t)=S(t)(\varphi(0)+g(\varphi))-g(u_{t}^{(n)})-\int_{0}^{t}AS(t-s)g(u_{s}^{n})ds+\int_{0}^{t}S(t-s)f(u_{s}^{n-1})ds\\\
+\int_{0}^{t}S(t-s)\sigma(u_{s}^{n-1})dW(s),\ t\in[0,T]$ (21)

with $u^{(n)}(t)=\varphi(t)$ for $t\in[-h,0]$. It was shown in [27] that if
the estimates (11) and (12) hold, the iteration scheme (21) converges to the
unique mild solution of (1). Hence the proof of Theorem 2.6 is complete in the
case, when $\mathbb{E}\|\varphi\|^{p}_{C_{h}}<\infty.$ In general case, we may
define a sequence of cut-off functions in spirit of [5], Th.7.4.

$\varphi_{n}:=\begin{cases}\varphi,\ \text{ if }\|\varphi\|_{C_{h}}\leq n\\\
0,\ \text{ if }\|\varphi\|_{C_{h}}>n.\end{cases}$

Let $u_{n}$ be the solution of (1) with the initial condition $\varphi_{n}.$
By [5], Th.7.4, the sequence $u_{n}$ converges to the process $u(t)$, which
solves (1). This completes the proof of Theorem 2.6.

###  3.2 Proof of Theorem 2.8.

Using Lemma 3.3, we have

$\mathbb{E}\sup_{s\in[0,t]}\|u^{(n)}(s)\|^{p}\leq
K_{1}+C_{1}\mathbb{E}\|\varphi\|^{p}_{C_{h}}+K\int_{0}^{t}\mathbb{E}\sup_{\tau\in[0,s]}\|u^{(n-1)}(s)\|^{p}d\tau.$

An analogous estimate holds for the solution $u(t)$ as well. Hence, by
Gronwall inequality we have

$\mathbb{E}\sup_{s\in[0,t]}\|u(s)\|^{p}\leq(K_{1}+C_{1}\mathbb{E}\|\varphi\|^{p}_{C_{h}})e^{Kt}.$
(22)

Next, for two solutions $u(t)$ and $v(t)$ with random initial conditions
$\varphi$ and $\psi$ respectively, using Lemma 3.3 we have

$\displaystyle\mathbb{E}\sup_{s\in[0,t]}\|u(s)-v(s)\|^{p}\leq
L\mathbb{E}\|\varphi-\psi\|^{p}_{C_{h}}+B\int_{0}^{t}\mathbb{E}N(\sup_{\tau\in[0,s]}\|u_{\tau}-v_{\tau}\|^{p})ds$
$\displaystyle\leq
L\mathbb{E}\|\varphi-\psi\|^{p}_{C_{h}}+BtN(\|\varphi-\psi\|^{p}_{C_{h}})+B\int_{0}^{t}\mathbb{E}N(\sup_{\tau\in[0,s]}\|u(\tau)-v(\tau)\|^{p})ds,$

where we implicitly used the fact that
$\sup_{\tau\in[-h,s]}N(\|u(\tau)-v(\tau)\|^{p})$ is attained either for
$\tau\in[-h,0]$ or for $\tau\in[0,s]$, and hence

$\mathbb{E}N(\sup_{\tau\in[-h,s]}\|u(\tau)-v(\tau)\|^{p})=\max\\{\mathbb{E}N(\sup_{\tau\in[-h,0]}\|u(\tau)-v(\tau)\|^{p}),\mathbb{E}N(\sup_{\tau\in[0,s]}\|u(\tau)-v(\tau)\|^{p})\\}\\\
\leq\mathbb{E}N(\sup_{\tau\in[-h,0]}\|u(\tau)-v(\tau)\|^{p})+\mathbb{E}N(\sup_{\tau\in[0,s]}\|u(\tau)-v(\tau)\|^{p}).$

It follows from (22) and that if
$\mathbb{E}\|\varphi_{n}-\varphi\|^{p}_{C_{h}}\to 0,n\to\infty$, we have

$\limsup_{n\to\infty}\mathbb{E}\sup_{s\in[0,t]}\|u(s,\varphi_{n})-u(s,\varphi)\|^{p}\leq
B\int_{0}^{t}N(\limsup_{n\to\infty}\mathbb{E}\sup_{\tau\in[0,s]}\|u(\tau,\varphi_{n})-u(\tau,\varphi)\|^{p})ds.$

Thus, from the condition (H2)(c), we have

$\limsup_{n\to\infty}\mathbb{E}\sup_{s\in[0,t]}\|u(s,\varphi_{n})-u(s,\varphi)\|^{p}=0,$

which completes the proof of Theorem 2.8.

###  3.3 Proof of Theorem 2.13.

The idea of the proof will be based on Krylov-Bogoliubov compactness approach
[5]. Let $u(t)$ be the solution of (1), which satisfies the condition (7). Our
goal is to show that the family of distributions $\mathcal{L}(u_{t}),t\geq T$,
with $T\geq 2h$, is tight. To this end, we will show that for any
$\varepsilon>0$ there is a compact $K_{\varepsilon}\subset C_{h}$, such that

$P\\{u_{t}\in K_{\varepsilon}\\}>1-\varepsilon.$ (23)

We have

$\displaystyle
u_{T}=u(T+\theta)=S(T+\theta)(\varphi(0)+g(\varphi))-g(u_{T+\theta})-\int_{0}^{T+\theta}AS(T+\theta-s)g(u_{s})ds$
(24)
$\displaystyle+\int_{0}^{T+\theta}S(T+\theta-s)f(u_{s})ds+\int_{0}^{T+\theta}S(T+\theta-s)\sigma(u_{s})dW(s).$

######  Lemma 3.4.

Let $\varphi(t,\omega):C_{h}\times\Omega\to C_{h}$ be a random process, such
that for any $\gamma>0$

$\lim_{\sigma\to
0}P\\{\omega:\sup_{t,s\in[-h,0],|t-s|<\sigma}\|\varphi(t,\omega)-\varphi(s,\omega)\|_{C_{h}}>\gamma\\}=0.$
(25)

Then for any $\varepsilon>0$ there exists $B_{\varepsilon}\in\Omega$ such
that:

  1. 1.

$P\\{B_{\varepsilon}\\}>1-\varepsilon$;

  2. 2.

the family of functions $\\{\varphi(t,\omega),\omega\in B_{\varepsilon}\\}$ is
equicontinuous on $[-h,0]$.

###### Proof.

It follows from (25) that for any $\varepsilon>0$ there exists
$\\{\delta_{k},k\geq 1\\}$ such that $\delta_{k}\to 0$ and

$P\\{\omega:\sup_{t,s\in[-h,0],|t-s|<\sigma}\|\varphi(t,\omega)-\varphi(s,\omega)\|_{C_{h}}>\frac{1}{k}\\}\leq\frac{\varepsilon}{2^{k}}.$
(26)

Denote

$A_{k}:=\\{\omega:\sup_{t,s\in[-h,0],|t-s|<\sigma}\|\varphi(t,\omega)-\varphi(s,\omega)\|_{C_{h}}>\frac{1}{k}\\}.$

It follows from (26) that

$P\\{\cup_{k\geq 1}A_{k}\\}\leq\sum_{k=1}^{\infty}P(A_{k})\leq\varepsilon.$

Set

$B:=\overline{\cup_{k\geq 1}A_{k}}=\cap_{k\geq 1}\overline{A_{k}}=\cap_{k\geq
1}\\{\omega:\sup_{t,s\in[-h,0],|t-s|<\sigma}\|\varphi(t,\omega)-\varphi(s,\omega)\|_{C_{h}}\leq\frac{1}{k}\\}.$

Then $P\\{B\\}>1-\varepsilon$, and if $\omega\in B$, then
$\omega\in\overline{A_{k}}$, and for any $\gamma>0$ choose
$k>\frac{1}{\gamma}$, which would imply
$\|\varphi(t,\omega)-\varphi(s,\omega)\|_{C_{h}}\leq\frac{1}{k}<\gamma$ for
$|t-s|<\delta_{k}$. This completes the proof of Lemma 3.4. ∎

######  Lemma 3.5.

Let $D$ be a set of functions $\varphi_{t}\in C([-h,0],C_{h})$ with the
following properties:

  1. 1.

there is $R>0$ such that

$\|\varphi_{0}\|_{C_{h}}\leq R\text{ for any }\varphi_{t}\in D.$

  2. 2.

the set of functions $D$ is equicontinous on $[-h,0]$ in $t$.

Then the set $\\{g(\varphi_{t}),\varphi\in D\\}$ is compact in $C_{h}$.

###### Proof.

We will apply the infinite dimensional version of Arzela-Ascoli Theorem. To
this end, we need to establish

  * [i]

For fixed $t\in[-h,0]$, $\\{g(\varphi_{t}),\varphi\in D\\}$ is compact in $H$;

  * [ii]

for any $\varepsilon>0$, there is $\delta>0$, such that

$\|g(\varphi_{t})-g(\varphi_{s})\|<\varepsilon,\text{ if }|t-s|<\delta,\ \
t,s\in[-h,0].$

It follows from the definition of $D$ that there is $C>0$ such that

$\sup_{t\in[-h,0]}\|\varphi_{t}\|_{C_{h}}\leq C\text{ for any }\varphi_{t}\in
D.$

Using (H3) we have

$\|g(\varphi_{t})\|_{\alpha}\leq M_{g}\|\varphi_{t}\|_{C_{h}}+\|g(0)\|\leq
M_{g}C+\|g(0)\|.$

This way Lemma 2.1 implies that the set $D$ is bounded in $H_{\alpha}$ and
therefore is compact in $H$. Furthermore, for all $\varepsilon>0$ we have

$\|g(\varphi_{t})-g(\varphi_{s})\|\leq\|g(\varphi_{t})-g(\varphi_{s})\|_{\alpha}\leq
M_{g}\|\varphi_{t}-\varphi_{s}\|_{C_{h}},$

which completes the proof of Lemma. ∎

######  Lemma 3.6.

The operator

$[B\varphi](\theta):=\int_{0}^{T+\theta}AS(T+\theta-s)\varphi(s)ds$

is a compact operator from $C([0,T],H_{\alpha})$ into $C_{h}$.

###### Proof.

We will also use the infinite dimensional analog of Arzela-Ascoli Theorem. For
fixed $\theta\in[-h,0]$ and $\varepsilon>0$ denote

$[B^{\varepsilon}\varphi](\theta):=\int_{0}^{T+\theta-\varepsilon}AS(T+\theta-s)\varphi(s)ds=\int_{0}^{T+\theta-\varepsilon}AS(\varepsilon)S(T+\theta-\varepsilon-s)\varphi(s)ds,$

Since $A$ is a generator of analytic semigroup, we have
$S(T+\theta-\varepsilon-s)\varphi(s)\in D(A)$, then

$\int_{0}^{T+\theta-\varepsilon}AS(\varepsilon)S(T+\theta-\varepsilon-s)\varphi(s)ds=S(\varepsilon)\int_{0}^{T+\theta-\varepsilon}AS(T+\theta-\varepsilon-s)\varphi(s)ds.$

But

$\left\|\int_{0}^{T+\theta-\varepsilon}AS(T+\theta-\varepsilon-s)\varphi(s)ds\right\|\leq\int_{0}^{T+\theta-\varepsilon}C_{1-\alpha}(T+\theta-\varepsilon-s)^{\alpha-1}ds\sup_{s\in[0,T]}\|\varphi(s)\|_{\alpha}\\\
\leq\frac{C_{1-\alpha}T^{\alpha}}{\alpha}\sup_{s\in[0,T]}\|\varphi(s)\|_{\alpha}.$

Since $S(\varepsilon)$ is compact, we conclude that the operators
$B^{\varepsilon}$ are compact. We have

$\|B^{\varepsilon}\varphi-B\varphi\|\leq\int_{T+\theta-\varepsilon}^{T+\theta}C_{1-\alpha}(T+\theta-s)^{\alpha-1}ds\sup_{s\in[0,T]}\|\varphi(s)\|_{\alpha}\to
0,\ \varepsilon\to 0.$

Consequently $B^{\varepsilon}$ converges to $B$ in the operator norm, thus $B$
is compact and the set

$\\{B\varphi(\theta):\sup_{t\in[0,T]}\|\varphi(t)\|_{\alpha}\leq 1\\}$

is compact in $H$. Fix $\theta$ and $r$, such that
$-h\leq\theta\leq\theta+r\leq 0$. Then for $\varphi$ satisfying
$\sup_{t\in[0,T]}\|\varphi(t)\|_{\alpha}\leq 1$ we have

$\|(B\varphi(\theta+r)-B\varphi(\theta)\|\leq\left\|\int_{0}^{T+\theta}(AS(T+\theta+r-s)-AS(T+\theta-s))\varphi(s)ds\right\|\\\
+\left\|\int_{T+\theta}^{T+\theta+r}AS(T+\theta+r-s)\varphi(s)ds\right\|:=J_{1}+J_{2}.$

In view of Lemma 2.2, the estimate for $J_{2}$ is straightforward:

$J_{2}\leq\int_{T+\theta}^{T+\theta+r}C_{1-\alpha}(T+\theta+r-s)^{\alpha-1}ds\sup_{t\in[0,T]}\|\varphi(t)\|_{\alpha}\to
0,r\to 0.$

In order to estimate $J_{1}$, we proceed as follows:

$J_{1}\leq\int_{0}^{T+\theta}\left\|(AS(s+r)-AS(s))\varphi(T+\theta-s)\right\|ds.$

Since $S(s)\varphi(T+\theta-s)\in D(A)$, using the continuity property of
semigroup, we have

$(AS(s+r)-AS(s))\varphi(T+\theta-s)=\\\
S(r)AS(s)\varphi(T+\theta-s)-AS(s)\varphi(T+\theta-s)\to 0,r\to 0.$

Furthermore,

$\|AS(s+r)\varphi(s)\|\leq\frac{C_{1-\alpha}}{s^{1-\alpha}}\sup_{t\in[0,T]}\|\varphi(t)\|_{\alpha}.$

Hence, by Dominated Convergence Theorem, $J_{1}\to 0$ as $r\to 0$, which
competes the proof of Lemma 3.6. ∎

The following result was shown in [21], Lemma 3.3:

######  Lemma 3.7.

For any $p>2$ and $\beta\in(\frac{1}{p},1]$ the operator

$(G_{\beta}\varphi)(\theta)=\int_{0}^{T+\theta}(T+\theta-s)^{\beta-1}S(T+\theta-s)\varphi(s)ds$

is a compact operator $L^{p}([0,T],H)$ to $C_{h}$.

For $\varphi\in C_{h}$, $\psi,\eta\in L^{p}([0,T],H)$, and $\xi\in
C([0,T],H)$, introduce

$\mu[\varphi,\xi,\psi,\eta]:=S(T+\theta)(\varphi(0)+g(\varphi))-g(u_{t})+\int_{0}^{T+\theta}AS(T+\theta-s)\xi(s)ds+(G_{1}\psi)(\theta)+(G_{\beta}\eta)(\theta),$

and denote

$K(r):=\\{\mu[\varphi,\xi,\psi,\eta]\in C_{h}:\|\varphi\|_{C_{h}}\leq
r,\sup_{t\in[0,T]}\|\xi\|\leq r,\|\psi\|_{L^{p}([0,T],H)}\leq
r,\|\eta\|_{L^{p}([0,T],H)}\leq r\\}.$

It follows from Lemma 3.6 and Lemma 3.7 that $K(r)$ is compact in $C_{h}$.
Let $u(t,\varphi)$ be a solution of (1) with the initial condition $\varphi\in
C_{h}$.

######  Lemma 3.8.

Introduce

$z(\theta):=S(T+\theta)(\varphi(0)+g(\varphi))+\int_{0}^{T+\theta}AS(T+\theta-s)g(u_{s}(\varphi))ds+\\\
\int_{0}^{T+\theta}S(T+\theta-s)f(u_{s}(\varphi))ds+\int_{0}^{T+\theta}S(T+\theta-s)\sigma(u_{s}(\varphi))dW(s).$

Assume the conditions of Theorem 2.6 hold. Then there is $C>0$ such that for
any $r>0$ and for any $\varphi\in C_{h}$ with $\|\varphi\|_{C_{r}}\leq r$, we
have

$P\\{z(\theta)\in K(r)\\}\geq 1-Cr^{-p}(1+\|\varphi\|^{p}_{C_{r}}).$

###### Proof.

Fix $p>2$ and $\beta\in\left(\frac{1}{p},\frac{1}{2}\right)$. Using the
factorization formula [4], Th.5.2.5, we have

$z(\theta)=S(T+\theta)(\varphi(0)+g(\varphi))+\int_{0}^{T+\theta}AS(T+\theta-s)g(u_{s}(\varphi))ds+\\\
(G_{1}f(u_{s}))(\theta)+\frac{\sin(\beta\pi)}{\pi}(G_{\beta}Y(s))(\theta),$

where

$Y(s)=\int_{0}^{s}(s-\tau)^{-\beta}S(s-\tau)\sigma(u_{\tau})dW(\tau).$

Lemma 7.2 [26] yields

$\mathbb{E}\int_{0}^{T}\|Y(s)\|^{p}ds=\mathbb{E}\int_{0}^{T}\left\|\int_{0}^{s}(s-\tau)^{-\beta}S(s-\tau)\sigma(u_{\tau})dW(\tau)\right\|^{p}\\\
\leq
C_{0}^{p}\mathbb{E}\int_{0}^{T}\left(\int_{0}^{s}(s-\tau)^{-2\beta}\|\sigma(u_{\tau}\|^{2}_{\mathcal{L}_{2}^{0}}d\tau\right)^{p/2}ds.$

Using Hausdorff-Young inequality and the condition (H2)(i), we get

$\mathbb{E}\int_{0}^{T}\|Y(s)\|^{p}ds\leq
C_{21}\left(\int_{0}^{T}t^{-2\beta}dt\right)^{\frac{p}{2}}\left(1+\mathbb{E}\sup_{t\in[0,T]}\|\varphi\|^{p}_{C_{h}}\right).$

Hence using (22) and Lemma 3.1 we have

$\mathbb{E}\int_{0}^{T}\|Y(s)\|^{p}ds\leq
C_{22}\left(1+\|\varphi\|^{p}_{C_{h}}\right).$

Similarly,

$\mathbb{E}\int_{0}^{T}\|f(u_{s})\|^{p}ds\leq
C_{23}\left(1+\mathbb{E}\|\varphi\|^{p}_{C_{h}}\right).$

Furthermore, using the condition (H3), (22) and Lemma 3.1 we have

$\mathbb{E}\sup_{s\in[0,T]}\|g(u_{s}(\varphi))\|^{p}_{\alpha}\leq
C_{24}\left(1+\|\varphi\|^{p}_{C_{h}}\right).$

This way

$P\\{z(\theta)\notin K(r)\\}\leq
P\\{\|f(u_{s}(\varphi))\|_{L^{p}([0,T],H)}>r\\}+P\\{\|Y(s)\|_{L^{p}([0,T],H)}>\frac{\pi
r}{\sin(\beta\pi)}\\}\\\
+P\\{\sup_{s\in[0,\tau]}\|g(u_{s}(\varphi))\|_{\alpha}>r\\}\leq
r^{-p}C\left(1+\mathbb{E}\|\varphi\|^{p}_{C_{h}}\right),$

which completes the proof of Lemma 3.8. ∎

######  Lemma 3.9.

The random process $u_{T+\theta}$ satisfies (25) if
$\mathbb{E}\|u_{0}\|^{p}_{C_{h}}<\infty$.

###### Proof.

For any $\gamma>0$ Chebyshev’s inequality yields

$P\\{\sup_{\theta_{1},\theta_{2}\in[-h,0]}\|u_{T+\theta_{1}}-u_{T+\theta_{2}}\|_{C_{h}}>\gamma\\}\leq\frac{1}{\gamma^{p}}\mathbb{E}\sup_{\theta_{1},\theta_{2}\in[-h,0]}\|u_{T+\theta_{1}}-u_{T+\theta_{2}}\|_{C_{h}}^{p}.$

Let $\mu:=\theta_{2}-\theta_{1}$. Then

$\sup_{\theta\in[-h,0]}\|u_{T+\theta+\mu}-u_{T+\theta}\|^{p}_{C_{h}}\leq\sup_{\theta\in[-h,0]}\sup_{s\in[T-h,T]}\|u(s+\theta+\mu)-u(s+\theta)\|^{p}\leq\sup_{s\in[0,T]}\|u(s+\mu)-u(s)\|^{p},$

since $T>2h.$ So $u(s,\omega)$ is continuous in $H$ with probability 1 on
$[0,T]$. Therefore $u(s)$ is uniformly continuous on $[0,T]$, so

$\sup_{s\in[0,T]}\|u(s+\mu)-u(s)\|^{p}\to 0\text{ a.s. as }\mu\to 0.$

Since

$\mathbb{E}\sup_{s\in[0,T]}\|u(s+\mu)\|^{p}\leq\mathbb{E}\sup_{s\in[0,T+1]}\|u(s)\|^{p}<\infty,$

by Dominated Convergence Theorem, we obtain

$\mathbb{E}\sup_{s\in[0,T]}\|u(s+\mu)-u(s)\|^{p}\to 0\text{ as }\mu\to 0,$

hence

$\mathbb{E}\sup_{\theta_{1},\theta_{2}\in[-h,0],|\theta_{1}-\theta_{2}|<\delta}\|u_{T+\theta_{1}}-u_{T+\theta_{2}}\|_{C_{h}}^{p}\to
0,\text{ as }\delta\to 0,$

which completes the proof of the Lemma. ∎

We now return to the proof of Theorem 2.13. For any $\varepsilon>0$ we will
show that there is a compact set $K_{\varepsilon}$ in $C_{h}$, which satisfies
the condition (23). Using Lemma 3.9, $u_{T+\theta}$ satisfies the condition
(25), hence there exists a set $B_{\varepsilon/3}$, such that
$P\\{B_{\varepsilon}\\}>1-\frac{\varepsilon}{3},$ and such that the set of
functions

$\\{u_{T+\theta}(\omega),\omega\in B_{\varepsilon/3}\\}$

is equicontinuous on $[-h,0]$. Next, choose $R_{1}>0$ such that

$P\\{\|u_{T}\|_{C_{h}}>R_{1}\\}\leq\frac{\varepsilon}{3}.$

Due to the condition (7), such choice is always possible. Therefore by Lemma
3.5 there is a compact set $K_{\frac{2\varepsilon}{3}}\in C_{h}$ such that

$P\\{g(u_{T+\theta})(\omega)\in K_{\frac{2\varepsilon}{3}}\\}\geq
1-\frac{2\varepsilon}{3}.$ (27)

Using the Markov property for $t>T$, for any Borel $K\subset C_{h}$ we have

$P\\{u_{t}\in K\\}=\mathbb{E}(P(u_{t}\in
K|\mathcal{F}_{t-1}))=\mathbb{E}(P(0,u_{t-T},T,K))\geq\mathbb{E}(P(0,u_{t-T},T,K)\chi_{\\{\|u_{t-T}\|^{p}_{C_{h}}\leq
R_{1}\\}}).$ (28)

Now, using Lemma 3.8, we can choose $r>R_{1}$ such that

$P\\{z(\theta)\in K(r)\\}\geq 1-\frac{\varepsilon}{3}.$ (29)

It follows from (24), (27) and (29) there is a compact $K_{\varepsilon}\subset
C_{h}$ such that

$P\\{u_{T}\in K(\varepsilon)\\}\geq 1-\varepsilon.$

Setting $K=K_{\varepsilon}$ in (28) we have

$P\\{u_{t}\in
K(\varepsilon)\\}\geq(1-\varepsilon)P\\{\|u_{t-T}\|^{p}_{C_{h}}\leq R_{1}\\}.$

In view of the condition (7), the proof of Theorem 2.13 is complete.

##  4 Applications.

In this section we use the standard elliptic operator $A$ to illustrate our
main result. We also provide explicit examples of non-Lipschitz nonlinearities
$f$ and $\sigma$, for which the main results are applicable.

Let $D$ be a bounded domain in $\mathbb{R}^{d}$ with $\partial D$ satisfying
the Lyapunov condition, $H=L^{2}(D)$, and

$Au=\sum_{i,j=1}^{d}(a_{i,j}(x)u_{x_{i}})_{x_{j}}={\rm{div}}(a(x)\nabla u).$

Here $a_{i,j}$ are Holder continuous coefficients with the Holder exponent
$\beta\in(0,1)$, symmetric, bounded, and satisfies the elliptic condition

$\sum_{i,j=1}^{d}a_{i,j}\eta_{i}\eta_{j}\geq
C_{0}\|\eta\|^{2},\eta\in\mathbb{R}^{d},$

for some $C_{0}>0.$ Let $e_{n}(x)$ be orthonormal basis in $H$, such that
$e_{n}\in L^{\infty}(D)$, and $\sup_{n}\|e_{n}\|_{L^{\infty}(D)}<\infty.$
Introduce the covariance operator $Q\in\mathcal{L}(H)$ such that $Q$ is non-
negative, $Tr(Q)<\infty$ and $Qe_{n}=\lambda_{n}e_{n}$. This enables us to
define

$W(t):=\sum_{i=1}^{\infty}\sqrt{\lambda_{i}}\beta_{i}(t)e_{i}(x),t\geq 0,$

which is a $Q$-Wiener process for $t\geq 0$ with values in $L^{2}(D)$. Denote
$U:=Q^{\frac{1}{2}}(L^{2}(D))$. It follows from [19], Lemma 2.2, that $U\in
L^{\infty}(D)$. Following [19], introduce the multiplication operator
$\Psi:U\to H$ as follows. For fixed $\varphi\in L^{2}(D)$, let
$\Psi(\varphi)=\varphi\cdot\psi$ for $\psi\in U$. Since $\varphi\in L^{2}(D)$
and $\psi\in L^{\infty}(D)$, the operator $\Psi$ is well defined, and hence
$\Psi\circ Q^{1/2}:L^{2}(D)\to L^{2}(D)$ defines a Hilbert-Schmidt operator,
with

$\|\Psi\circ Q^{1/2}\|^{2}_{\mathcal{L}_{2}^{0}}\leq
Tr(Q)\sup_{n}\|e_{n}\|^{2}_{\infty}\|\varphi\|^{2}_{L^{2}(D)}.$

Our main object of interest in this section is the following delay equation:

$d[u(t,x)+\int_{D}b(x,u(t-h,y),y)dy]=[{\rm div}(a(x)\nabla
u(t,x))+f(u(t-h,x))]dt+\sigma(u(t-h,x))dW(t)$ (30)

for $t>0$, with $u(t,x)=\varphi(t,x)$ for $t\in[-h,0]$ and $u(t,x)=0$ for
$x\in\partial D,t\geq 0.$ Here
$b(x,z,y):\mathbb{R}^{d}\times\mathbb{R}^{1}\times\mathbb{R}^{d}\to\mathbb{R}$,
$f:\mathbb{R}\to\mathbb{R}$ and $\sigma:\mathbb{R}\to\mathbb{R}.$ Introduce
the mapping

$g(\varphi)(x):=\int_{D}b(x,\varphi(-h),y)dy$

as a map from $C_{h}$ to $L^{2}(D)$. Then the problem (30) can be viewed in
the abstract form (1), with $D(A)=H^{2}(D)\cap H_{0}^{1}(D)$.

######  Corollary 4.1.

  * [i]

Suppose $f$ and $\sigma$ satisfy (H2) (as Nemytskii maps), and the initial
condition $\varphi$ satisfies (H4). In addition, assume the function $b$ is
continuous with respect to all of its variables, and there exist constants
$A>0$ and $M_{g}>0$ such that

$|b(x,0,y)|+|\nabla_{x}b(x,z,y)|\leq A$

and

$|b(x,z_{1},y)-b(x,z_{2},y)|+|\nabla_{x}b(x,z_{1},y)-\nabla_{x}b(x,z_{2},y)|\leq
M_{g}|z_{1}-z_{2}|$

for all $x,y\in D$, $z_{1},z_{2}\in\mathbb{R}$ and

$2M_{g}^{2}\ {\rm meas}^{2}(D)<1.$ (31)

Then (30) has a unique mild solution.

  * [ii]

If, in addition, the functions $f$, $\sigma$ and $b$ are bounded, there is an
invariant measure for (30).

###### Proof.

Let us verify that the conditions of Theorem 2.6 hold. Since $(-A)$ is a
sectorial self-adjoint operator, it follows, e.g. from [7], p.335 that
$\sigma(-A)>\sigma_{0}>0$, and $(-A)^{-1}$ is compact. Therefore the
assumption (H1) holds. Furthermore, following [9] introduce the interpolation
space $D_{A}(\frac{1}{2},2)=H_{0}^{1}$. Therefore by Proposition A.17 [14]
$D_{A}(\frac{1}{2},2)$ is isomorphic to $D((-A)^{\frac{1}{2}})$, and

$\|g(\varphi)\|^{2}_{\frac{1}{2}}=\|g(\varphi)\|^{2}_{H_{0}^{1}}=\int_{D}|(g(\varphi)(x)|^{2}dx+\int_{D}|\nabla
g(\varphi)(x)|^{2}dx.$

But

$\int_{D}|g(\varphi)(x)|^{2}dx=\int_{D}dx\left(\int_{D}b(x,\varphi(-h,y),y)dy\right)^{2}\leq
L^{2}{\rm meas}^{2}(D)\left(\int_{D}|\varphi(-h,y)|^{2}dy+1\right)<\infty.$

By Lebesgue Theorem on differentiation of integrals,

$\nabla g(\varphi)(x)=\int_{D}\nabla_{x}b(x,\varphi(-h,y),y)dy,$

hence

$\int_{D}|\nabla g(\varphi)(x)|^{2}dx<\infty.$

Finally,

$\displaystyle\|g(\varphi_{1})-g(\varphi_{2})\|^{2}_{\frac{1}{2}}=\int_{D}|g(\varphi_{1})(x)-g(\varphi_{2})(x)|^{2}dx+\int_{D}|\nabla_{x}g(\varphi_{1})(x)-\nabla_{x}g(\varphi_{2})(x)|^{2}dx$
$\displaystyle\leq\int_{D}dx\left(\int_{D}|b(x,\varphi_{1}(-h,y),y)-b(x,\varphi_{2}(-h,y),y)|dy\right)^{2}$
$\displaystyle+\int_{D}dx\left(\int_{D}|\nabla_{x}b(x,\varphi_{1}(-h,y),y)-\nabla_{x}b(x,\varphi_{2}(-h,y),y)|dy\right)^{2}$
$\displaystyle\leq 2{\rm
meas}^{2}(D)M_{g}^{2}\sup_{\theta\in[-h,0]}\int_{D}|\varphi_{1}(\theta,y)-\varphi_{2}(\theta,y)|^{2}dy=2{\rm
meas}^{2}(D)M_{g}^{2}\|\varphi_{1}-\varphi_{2}\|^{2}_{C_{h}},$

which implies that the condition (H3) holds. Thus [i] follows by Theorem 2.6.

In order to derive [ii], it remains to verify the condition (7) and apply
Theorem 2.13. To this end, the spectral condition on $A$ implies the
exponential contraction property of the semigroup

$\|S(t)\|\leq Ke^{-\sigma_{0}t},t\geq 0.$

Assume

$|f|^{2}+|\sigma|^{2}+|b|^{2}+|\nabla_{x}b|^{2}\leq C.$

Then by Lemma 3.1

$\sup_{t\geq
0}\mathbb{E}\|u_{t}\|^{2}_{C_{h}}\leq\mathbb{E}\|\varphi\|^{2}_{C_{h}}+\sup_{t\geq
0}\mathbb{E}\|u(t)\|^{2}_{L^{2}(D)},$

which implies

$\displaystyle\mathbb{E}\|u(t)\|^{2}\leq
5\mathbb{E}\|S(t)(\varphi(0)+g(\varphi))\|^{2}+5\mathbb{E}\|g(u_{t})\|^{2}+5\mathbb{E}\left\|\int_{0}^{t}AS(t-s)g(u_{s})ds\right\|^{2}$
$\displaystyle+5\mathbb{E}\left\|\int_{0}^{t}S(t-s)f(u_{s})ds\right\|^{2}+5\mathbb{E}\left\|\int_{0}^{t}S(t-s)\sigma(u_{s})dW(s)\right\|^{2}\leq
5K^{2}e^{-2\sigma_{0}t}\mathbb{E}\|\varphi(0)+g(\varphi)\|^{2}$ $\displaystyle
C^{2}{\rm
meas}^{2}(D)+\mathbb{E}\left(\int_{0}^{t}\left\|A^{1/2}S(t-s)\|\|g(u_{s})\|_{1/2}\right\|ds\right)^{2}$
$\displaystyle
K^{2}\int_{0}^{t}e^{-\sigma_{0}(t-s)}Cds+K^{2}\sum_{k=1}^{\infty}\lambda_{k}\sup_{n\geq
1}\|e_{n}\|^{2}_{\infty}\int_{0}^{t}e^{-\sigma_{0}(t-s)}ds.$

Since $\|g(u_{s})\|_{1/2}<\infty$ and $\|A^{1/2}S(t-s)\|\leq
C_{1/2}(t-s)^{-\frac{1}{2}}e^{-\sigma_{0}(t-s)}$, we have $\sup_{t\geq
0}\mathbb{E}\|u(t)\|^{2}<\infty$. Therefore, the condition (7) holds by
Chebyshev’s inequality. ∎

######  Example 4.2.

Suppose $meas(D)=1$. For $x\geq 0$ let

$N(x):=\begin{cases}0,\ x=0;\\\ -x\ln(x),\ 0<x\leq e^{-2};\\\ x+e^{-2},\
x>e^{-2}.\end{cases}$ (32)

In view of Remark 2.4, this choice of $N$ satisfies the conditions (a) and (c)
in (H2). For $p>2$ let

$f(x):=\begin{cases}0,\ x=0;\\\ |x|p^{1/p}\left|\ln|x|\right|^{1/p},\
0<|x|\leq e^{-2};\\\ \text{ any bounded Lipschitz function with }L\leq
1,\,|x|>e^{-2}.\end{cases}$ (33)

We first claim that for all $x,y\in[0,e^{-2}]$ we have

$|f(x)-f(y)|^{p}\leq N(|x-y|^{p})$ (34)

Indeed, since $f$ is concave on $[0,e^{-2}]$ and $f(0)=0$, it is also sub-
additive on this interval, i.e. $\forall a,b:f(a)+f(b)\geq f(a+b)$. The latter
inequality means that if $0<x<y<e^{-2}$, we have

$f(x)+f(y-x)\geq f(y),$

or

$(f(y)-f(x))^{p}\leq f^{p}(y-x)=N\left((y-x)^{p}\right).$

Using the same reasoning if $y<x$, the bound (34) follows. Since $f$ is even
on $[-e^{-2},e^{-2}]$ and Lipschitz outside of this interval, it is easy to
see that in fact (34) holds for all $x,y\in\mathbb{R}$. Note that since
$meas(D)=1$, the bound (34) is equivalent to the inequality

$\|f(\varphi_{1})-f(\varphi_{2})\|^{p}\leq N(\|\varphi_{1}-\varphi_{2}\|^{p})$
(35)

which, in turn, yields (H2) (b) for two constant initial conditions
$\varphi_{1}(t,x)\equiv\varphi_{1}$ and $\varphi_{2}(t,x)\equiv\varphi_{2}$,
$x\in D$, $t\in[-h,0]$. In order to obtain the estimate (H2) (b) for non-
constant $\varphi_{1}$ and $\varphi_{2}$, we equi-partition $D$ into
subdomains, approximate $\varphi_{1}$ and $\varphi_{2}$ with constants in each
subdomain, and use (34) in conjunction with the concavity of $N$ to get (35)
for piecewise constant functions. The desired result will then follow if we
take the limit as the number of partitions goes to infinity.

Thus, the equation (30) with a non-Lipschitz $f$, given by (33), is an example
of a well-posed neutral type stochastic functional-differential equation,
which admits the existence of invariant measure. The same conclusion will hold
if we choose $\sigma\equiv f$, given by (33). The only difference in this case
that we need to slightly modify the choice of $N$ in (32) by taking
$\tilde{N}:=\sum_{k=1}^{\infty}\lambda_{k}\sup_{n}\|e_{n}\|_{\infty}^{2}N$.

## Acknowledgement

The research of Oleksandr Misiats was supported by Simons Collaboration Grant
for Mathematicians No. 854856. The work of Oleksandr Stanzhytskyi was
partially supported by the National Research Foundation of Ukraine No.
F81/41743 and Ukrainian Government Scientific Research Grant No. 210BF38-01.

## References

  * [1] A.Anguraj and A. Vinodkumar.  Existence, uniqueness and stability of impulsive stochactic partial neutral functional differential equations with infinite delays.  J. Appl. Math. Informatics, 28(3-4).
  * [2] B. Boufoussi, S. Hajji, and E. Lakhel.  Time-dependent neutral stochastic functional differential equations driven by a fractional Brownian motion.  Commun. Stoch. Anal., 10(1):Article 1, 1–12, 2016.
  * [3] Carlo Cattaneo.  Sulla conduzione del calore.  Atti Sem. Mat. Fis. Univ. Modena, 3:83–101, 1949.
  * [4] G. Da Prato and J. Zabczyk.  Ergodicity for infinite-dimensional systems, volume 229 of London Mathematical Society Lecture Note Series.  Cambridge University Press, Cambridge, 1996.
  * [5] Giuseppe Da Prato and Jerzy Zabczyk.  Stochastic equations in infinite dimensions, volume 44 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 1992.
  * [6] D. A. Dawson.  Stochastic evolution equations.  Math. Biosci., 15:287–316, 1972.
  * [7] Lawrence C. Evans.  Partial differential equations, volume 19 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, second edition, 2010.
  * [8] T. E. Govindan.  Mild solutions of neutral stochastic partial functional differential equations.  Int. J. Stoch. Anal., pages Art. ID 186206, 13, 2011.
  * [9] P. Grisvald.  Commutativite de deux foncteurs d’interpolation et applications.  J. Math. Pures. Appl., 45:193–290, 1966.
  * [10] Morton E. Gurtin and A. C. Pipkin.  A general theory of heat conduction with finite wave speeds.  Arch. Rational Mech. Anal., 31(2):113–126, 1968.
  * [11] Daniel Henry.  Geometric theory of semilinear parabolic equations, volume 840 of Lecture Notes in Mathematics.  Springer-Verlag, Berlin-New York, 1981.
  * [12] Matthias Hieber, Oleksandr Misiats, and Oleksandr Stanzhytskyi.  On the bidomain equations driven by stochastic forces.  Discrete Contin. Dyn. Syst., 40(11):6159–6177, 2020.
  * [13] Feng Jiang and Yi Shen.  A note on the existence and uniqueness of mild solutions to neutral stochastic partial functional differential equations with non-Lipschitz coefficients.  Comput. Math. Appl., 61(6):1590–1594, 2011.
  * [14] K. Kenzhebaev, O. Stanzhytskyi, and A. Tsukanova.  Existence and uniqueness results, the markovian property of solutions for a neutral delay stochastic reaction-diffusion equation in entire space.  Dynamic Systems and Applications, 28(1):19–46, 2019.
  * [15] Nicolas Kryloff and Nicolas Bogoliouboff.  La théorie générale de la mesure dans son application à l’étude des systèmes dynamiques de la mécanique non linéaire.  Ann. of Math. (2), 38(1):65–113, 1937.
  * [16] Qiang Li and Mei Wei.  Existence and asymptotic stability of periodic solutions for neutral evolution equations with delay.  Evolution Equations and Control Theory, 9(3).
  * [17] Kai Liu.  Stationary solutions of neutral stochastic partial differential equations with delays in the highest-order derivatives.  Discrete Contin. Dyn. Syst. Ser. B, 23(9):3915–3934, 2018.
  * [18] Jiaowan Luo.  Exponential stability for stochastic neutral partial functional differential equations.  J. Math. Anal. Appl., 355(1):414–425, 2009.
  * [19] Ralf Manthey and Thomas Zausinger.  Stochastic evolution equations in $L^{2\nu}_{\rho}$.  Stochastics Stochastics Rep., 66(1-2):37–85, 1999.
  * [20] James Clerk Maxwell.  On the dynamical theory of gases.  Philosophical transactions of the Royal Society of London, (157):49–88, 1867.
  * [21] Oleksandr Misiats, Vikotria Mogilova, and Oleksandr Stanzhytskyi.  Invariant measure for stochastic functional differential equations in hilbert spaces.  arXiv:2011.07034, 2020.
  * [22] Oleksandr Misiats, Oleksandr Stanzhytskyi, and Nung Kwan Yip.  Existence and uniqueness of invariant measures for stochastic reaction-diffusion equations in unbounded domains.  J. Theoret. Probab., 29(3):996–1026, 2016.
  * [23] Oleksandr Misiats, Oleksandr Stanzhytskyi, and Nung Kwan Yip.  Asymptotic analysis and homogenization of invariant measures.  Stoch. Dyn., 19(2):1950015, 27, 2019.
  * [24] Oleksandr Misiats, Oleksandr Stanzhytskyi, and Nung Kwan Yip.  Invariant measures for stochastic reaction-diffusion equations with weakly dissipative nonlinearities.  Stochastics, 92(8):1197–1222, 2020.
  * [25] A. Pazy.  Semigroups of linear operators and applications to partial differential equations, volume 44 of Applied Mathematical Sciences.  Springer-Verlag, New York, 1983.
  * [26] Dehao Ruan and Jiaowan Luo.  The existence, uniqueness, and controllability of neutral stochastic delay partial differential equations driven by standard Brownian motion and fractional Brownian motion.  Discrete Dyn. Nat. Soc., pages Art. ID 7502514, 11, 2018.
  * [27] A. M. Samoilenko, N. I. Mahmudov, and A. N. Stanzhitskii.  Existence, uniqueness, and controllability results for neutral FSDES in Hilbert spaces.  Dynam. Systems Appl., 17(1):53–70, 2008.
  * [28] O. Stanzhytskyi.  Investigation of invariant sets of Itô stochastic systems by means of Lyapunov functions.  Ukraïn. Mat. Zh., 53(2):282–285, 2001.
  * [29] O. Stanzhytskyi.  Investigation of the exponential dichotomy of Itô stochastic systems by means of quadratic forms.  Ukraïn. Mat. Zh., 53(11):1545–1555, 2001.
  * [30] Oleksandr Stanzhytskyi, Viktoria V. Mogilova, and Alisa O. Tsukanova.  On comparison results for neutral stochastic differential equations of reaction-diffusion type in $L_{2}(\mathbb{R}^{d})$.  In Modern mathematics and mechanics, Underst. Complex Syst., pages 351–395. Springer, Cham, 2019.
  * [31] E. F. Tsarkov.  Random Perturbations of Functional Differential Equations.  “Zinatne”, Riga, 1989.


# Study of a large array to detect ultra high energy tau-neutrino

S. ATIK YILMAZ<EMAIL_ADDRESS>Bolu Abant Izzet Baysal University,
Bolu 14280, TR  A.YILMAZ<EMAIL_ADDRESS>Giresun University, Giresun
28200, TR  M. IORI<EMAIL_ADDRESS>La Sapienza, University of
Rome, Rome 00185, IT  K. Y. OYULMAZ<EMAIL_ADDRESS>Bolu Abant Izzet
Baysal University, Bolu 14280, TR  H. DENIZLI<EMAIL_ADDRESS>Bolu Abant
Izzet Baysal University, Bolu 14280, TR

###### Abstract

The $PeV$ cosmogenic neutrino is still interesting argument. Since cosmogenic
neutrinos interact weakly with matter, the detection of their direction will
precisely point out the source in the space. In this paper, we show the
results of the simulation of tau lepton air showers induced by high energy
neutrinos detected by an array of stations designed to use the Earth Skimming
method improved by the “mountain chain screen” strategy. Both track time stamp
and position information of the stations on the array are used to reconstruct
the shower to estimate the direction and the number of events. The array
studied consists of 640 stations ($40\times 16$) spread over an area of
$0.6\,km^{2}$ starting from $1500\,m$ above the sea level (a.s.l.) on
$30\degree$ inclined plane of the mountain. When we extrapolate to 3 years and
10 $km^{2}$ we estimate 13 tau lepton events in energy interval of 10 PeV to
1000 PeV detected using the present upper limits of tau neutrino flux.

Keywords: Ultra High Energy Cosmic Rays, Tau-Neutrino, Ground Based Array,
Earth Skimming Strategy, CORSIKA

##  1 Introduction

The observation of Ultra High Energy (UHE) neutrinos generated in both distant
galactic and extra galactic sources such as pulsars or supernova remnants have
been a major challenge in astroparticle physics. The neutrinos with an energy
of $10^{17}$ eV or more, generated by interaction of Ultra High Energy Cosmic
Rays (UHECRs) with the Cosmic Microwave Background (CMB) are called cosmogenic
neutrinos [1]. Because of the low neutrino flux and the small detection
probability, a very large scale detector array (surface or volume) is required
to get a detectable rate. For this purpose, many large array experiments based
on different techniques have been constructed [2, 3, 4, 5, 6, 7, 8].

Recently the IceCube Collaboration [9] has reported 28 neutrino candidate
events in the energy range of $30$ to 1200 TeV. They also claimed the
detection of a high-energy neutrino event, IceCube-170922A, with an energy of
$\sim 290$ TeV in correlation with other telescopes showing a multimessenger
physics. The arrival direction of this neutrino was consistent with the
location of a known $\gamma$-ray blazar, TXS 0506+056, observed to be in a
flaring state [10].

The cosmogenic neutrino, such as tau neutrino ($\nu_{\tau}$), may interact
along the chord of the Earth or inside a mountain chain and creates a tau
lepton ($\tau$) by the charged-current (CC) interaction. When the $\tau$
emerges from the rock before it decays, it produces an extensive air shower
developing at very large angle near to the horizon [11, 12]. The observations
can be achieved using surface detectors, Cherenkov light telescopes,
fluorescence detectors or radio signals. Some of the proposed experiments in
the literature are Ashra-1 [13], MAGIC telescopes [14], Telescope Array (TA)
[15], and IceCube [16, 17].

The Ashra-1 detector aims to observe the Cherenkov lights from $\tau$ showers
with Earth Skimming method for PeV-EeV energy range [13]. An upper limit of
neutrino flux has been measured by MAGIC telescopes in the 1 PeV to 3 EeV
energy interval for $\nu_{\tau}$ induced showers arising from the ocean [14].
In Antartica, radio antennas are placed to the ground and point at the high
mountains to detect the skimming $\nu_{\tau}$s too [18]. All these proposed
experiments have no evidence of $\tau$ events except a possible $\tau$
candidate in IceCube [16, 17]. To detect the horizontal and upward
$\nu_{\tau}$ showers, we propose a surface detector array on an inclined plane
to improve the detection acceptance of $\tau$ shower using the time of flight
(TOF) and e/mu separation [19, 20]. In this study we also include a chain of
mountain in front of the array and evaluate the wideness of the valley.

This paper shows the results of a large array named as “TAUshoWER”, (TAUWER)
which has the geometrical advantage of using both the Earth Skimming method to
detect the $\tau$ showers produced below the horizon and a mountain chain in
front of the array to enlarge the acceptance. The array geometry is discussed
in Section 2 while the simulation of the decaying $\tau$ air showers, the
selection criteria for decay length from emerging point and the result of the
probability calculations for tau-neutrino interaction are given in Section 3.
Analytical approach on the estimation of the event number to be observed with
this array is presented in Section 3.1. Reconstruction of the showers is given
in Section 3.2 and Section 3.3. Section 3.4 contains a discussion of the
trigger. Finally, we summarize our results and conclude with number of
expected event/year.

##  2 Design of the Array

To detect the $\tau$ showers produced by the tau-neutrino interaction with
rock along its path in the Earth crust or a mountain in front of the array, we
propose a large surface array in which stations point below the horizon to
detect particles produced by large angle showers. The best way, as discussed
in Section 3, to detect $\nu_{\tau}$ is to use large amount of matter by
pointing down the horizon and a mountain chain in front of the array to
increase the interaction probability of $\nu_{\tau}$. Therefore, we propose to
install stations on an inclined mountain surface. In this paper, we consider a
grid of 640 stations ($40\times 16$) placed $30\,m$ apart and located on
$30\degree$ inclined plane of the mountain between $1500\,m$ and $2250\,m$
above the sea level and covering a surface of 0.6 km2 as shown in Figure 1.
Each station, named tower, is composed by two scintillator tiles ($40\times
20\times 1.5\,cm^{3}$) 160 cm apart and read by 3x3 mm2 silicon
photomultiplier (SiPM) placed on one side , which permit to select the
particle direction using the time of flight (TOF) method. If a layer of lead
2.5 cm thickness is installed in front of the rear tile it is possible
separate the muons from the electromagnetic particles [20]. The stations point
below the horizon (approximately $2.5\degree$) to measure upward moving
showers coming from $\nu_{\tau}$ interactions with the Earth crust and at the
some time can detect showers produced in the mountain chain. Shower direction
is selected by the TOF computed with the two scintillating plates. The layout
of the array studied in this paper covers a surface of 0.6 km2 which can be
increased according to the shape of the mountain.

Figure 1: Schematic drawing (scaled) of the array. The array consists of 640
stations placed on the $30\degree$ inclined plane (mountain surface) in the
matrix form of 16 rows and 40 columns which are separated by 30 m apart.

##  3 Evaluation of $\tau$ events escaped from the rock

The detection of high energetic ${\nu}_{\tau}$ which are coming nearly
horizontal to the Earth surface has a chance to create a $\tau$ via a charge-
current (CC) interaction or neutral-current (NC) interaction. It may survive
and emerge from the surface to initiate a shower in the atmosphere. In order
to take advantage of both valley and mountain to improve the detection
probability, each station in the array is placed on the mountain slope and
directed to the horizon. The shower detection probability depends on not only
$\tau$ initial energy but also traveling distance of the $\tau$ after the
$\nu_{\tau}$ interaction inside the Earth crust. Therefore, we have considered
a mountain in front with 60 km thickness which corresponds to maximum travel
distance for $\tau$ with an energy of 1000 PeV to escape from the rock before
it decays as shown in the sketch of the proposed experiment for this study
(Figure 2). The incoming $\nu_{\tau}$ travels a distance $(L-x)$ along the
Earth crust and interacts with the rock as depicted in the asterisk symbol,
$``*"$ and produce $\tau$ by CC interaction. The $\tau$ travels a distance x
before emerging from the rock. $L^{*}$ is the distance between the point where
$\tau$ initiates the shower and the center of the array. In this study, we
have required that $\tau$ shower is produced in atmosphere.

Figure 2: In this layout (not in scale) the array accepts $\tau$ shower
produced by charged current interaction skimmed neutrinos by the Earth crust
and a mountain in front of the array. $L-x$ is the interaction distance for
the ${\nu}_{\tau}$, the interaction point shown by an asterisk symbol $``*"$,
$x$ is the traveling distance for the ${\tau}$. $L^{*}$ is the distance of
$\tau$ decay point from the center of the array.

To evaluate the expected $\tau$ events detected by the array we first
generated events using the TAUOLA code [21] considering several decay modes of
the $\tau$ in the energy range of 10-1000 PeV. We have considered the decay
modes ($\pi^{-}$, $\pi^{-}\pi^{0}$, $\pi^{-}\pi^{+}\pi^{-}$,
$\pi^{-}\pi^{0}\pi^{0}$, $\pi^{-}\pi^{+}\pi^{-}\pi^{0}$) that covers 64 % of
the $\tau$ branching ratio to study the optimization and identification
performance of the expected $\nu_{\tau}$ induced showers by TAUWER array.

Consequently, CORSIKA (version 6.99) [22] has been modified and compiled to
simulate very inclined showers at an observation plane which is considered to
be $30\degree$. The produced air showers, induced by the decay products of
$\tau$, were initiated at 1500 m a.s.l. and developed up to the detector level
(2250 m a.s.l.) for 4 different distances starting from the $\tau$ decay point
to the center of the array, $L^{*}$ (3 km, 5 km, 7 km and 10 km), as shown in
Figure 2. 1000 showers have been simulated at different energy for this study.
QGSJETII [23] and GEISHA [24] are selected for high and low energy interaction
models, respectively. “CURVED EARTH” and “SLANT” options are activated to make
more realistic simulation for correctly include the atmospheric depth. “COAST”
option is selected for reading and converting CORSIKA binary files. Since the
computing time of the showers increases with the energy, the “thinning”
strategy has been also used in shower productions. The CORSIKA output file
provides position ($x$, $y$ and $z$), arrival time and momentum ($P_{x}$,
$P_{y}$, $P_{z}$, and $P_{tot}$) informations of each particle crossing on the
observation plane. Using these informations the trajectory is reconstructed
for each particle and counted when it passes trough the scintillating tiles.
In this proposed layout it is also important to identify the best distance of
the plane of the array from the $\tau$ decay point to intercept the maximum
shower evolution to make more efficient the shower detection and in the some
time evaluate the best valley wideness.

Figure 3: Averaged number of ${e^{\pm}}$ per shower (among 180 simulated
showers) detected on the array as a function of $L^{*}$, distance of $\tau$
decay point to the center of the array, for different energies ($10^{7}$,
$5x10^{7}$, $10^{8}$, $5x10^{8}$, $10^{9}$ GeV) at $\pi^{-}\pi^{0}$ decay
mode.

The determination of the optimum wideness for the valley is evaluated using
the particle density measured on the array plane for different decay paths of
the $\tau$ in air. Figure 3 shows the averaged number of the ${e^{\pm}}$ per
shower detected on the array plane for $\pi^{-}\pi^{0}$ modes at different
energies (10 - 1000 PeV) as a function of $L^{*}$, the distance of the $\tau$
decay point to the center of the array. The number of muons is about 5$\%$ of
the number of electrons and almost flat as a function of the distance of
$\tau$ decay point in the air. It shows that the number of $e^{\pm}$ from the
$\pi^{-}\pi^{0}$ decay mode and different $\tau$ energy on the detector plane
is maximum when the $\tau$ decay point is at a distance of 5-7 km from the
center of the array. A similar distribution is also obtained by $\pi^{-}$,
$\pi^{-}\pi^{+}\pi^{-}$, $\pi^{-}\pi^{0}\pi^{0}$ and
$\pi^{-}\pi^{+}\pi^{-}\pi^{0}$ decay modes. The decrement on the density of
$e^{\pm}$ at 10 km is due to the energy loss of $e^{\pm}$ in the atmosphere.
Since the density of the detected particles is higher at 5-7 km, 7 km decay
lengths from emerging point to the array center is used for further analysis
in this text. To obtain the expected number of neutrino $\tau$ detected by the
array we have to evaluate the probability of escaping of the $\tau$ from the
Earth crust, $P_{escp}$. The probability of $\nu_{\tau}$ interacts along the
$L-x$ path at $*$ point inside the rock is given as

$P_{int}=\int_{0}^{L-x}e^{-l/L_{int}}\frac{dl}{L_{int}}$ (1)

where $L_{int}$ is the neutrino interaction length, $dl$ is the short distance
the interaction happens. The neutrino cross section for charged current
interaction within 10$\%$ as follows according to CTEQ6.6M [25] and its energy
dependence varies as ${E_{\nu}}^{0.348}$ where $E_{\nu}$ is the $\nu_{\tau}$
energy. By this parameterization we have an interaction length of $\tau$ about
700 km for 10 PeV and 3000 km for 200 PeV. The neutrino NC cross section is
2.4 smaller than that of CC interaction which means the interaction length of
neutrino is about two times greater than that of charged current channel.
Since the NC interaction probability is no high before the CC interaction, the
neutrino has the same energy as it enters the Earth shell when it interacts as
CC. In this study we consider CC interactions of $\nu_{\tau}$ where the $\tau$
lepton is produced with 80% of the $\nu_{\tau}$ energy [26] as well as the NC
interaction.

On the other hand, the probability the $\tau$ escapes from the Earth shell is
given by the probability the neutrino interacts along its path times the
survival probabilities of $\tau$ lepton in the shell with an initial energy as
a function of the distance it travels, given as

$P_{escp}=P_{int}\times
exp\Bigg{[}-\displaystyle\int_{L-x}^{L}(ct_{0}\,(exp[-\beta_{0}/\beta_{1}+(\beta_{0}/\beta_{1}+\ln{(E_{i}/E_{0}}))]exp[-\beta_{1}\times\rho_{rock}\times
l]E_{0})/m_{\tau})^{-1}dl\Bigg{]}$

(2)

where $l$ is the distance the $\tau$ travelled, $c$ is the speed of light,
$t_{0}$ is the mean life time of $\tau$ lepton, $\beta_{0}$ is the energy loss
parameter from Bremsstrahlung, pair production and photonuclear interaction
with $E_{0}=10^{10}$ GeV, $E_{i}$ is the energy of $\tau$ lepton when it is
generated, $\rho_{rock}=2.65\times 10^{15}$ (g/km3) is the density of the
standard rock, and $m_{\tau}$ is the mass of $\tau$ lepton, $\beta_{1}$ is the
photonuclear coefficient as function of neutrino energy [27, 25, 28]. The
maximum distance the $\tau$ lepton travels at the energy of 1000 PeV is about
60 km. That means the neutrino interaction should happen in the Earth shell of
60 km thickness to produce a tau lepton which has a probability to escape from
the Earth crust. The survival probability ($P_{surv}$) of $\tau$, is
correspond to the exponential factor in Eq. 2, decreases while $\tau$ lepton
moves away from the interaction point in the Earth shell.

###  3.1 Estimation on the number of detected ${\nu}_{\tau}$ events

As discussed in Section 3, we evaluated the $\nu_{\tau}$ interaction length
and its probability $P_{int}$ based on the study given in Ref.[26], besides
that the approaches are used from Ref.[29] to calculate the $\tau$’s
$P_{surv}$ and then the escape probability, $P_{escp}$ after the $\tau$
travels through the Earth and the mountain chain $60$ km thickness at PeV
energies. Table 1 shows the $P_{escp}$ of the $\tau$ calculated for different
neutrino initial energies and escaping directions from the Earth crust using
the formula 2. The $\tau$ crossing the mountains are produced at $90\degree$
in $\theta$. No specific profile of the mountain has been studied.

Table 1: The escape probability of $\tau$ as function of $\theta$ and
$\nu_{\tau}$ energy. In the first row is reported the probability considering
a mountain screen 60 km thickness. These values where computed using a
neutrino cross section for charged currents according to CTEQ6.6M.

$\theta$(°) / Energy (PeV) 10 20 50 100 200 300 500 1000 Screen $6.850\times
10^{-4}$ $3.168\times 10^{-3}$ $9.430\times 10^{-3}$ $1.546\times 10^{-2}$
$2.220\times 10^{-2}$ $2.651\times 10^{-2}$ $3.245\times 10^{-2}$ $4.164\times
10^{-2}$ $91$ $3.114\times 10^{-6}$ $2.038\times 10^{-5}$ $9.056\times
10^{-5}$ $1.900\times 10^{-4}$ $3.276\times 10^{-4}$ $4.203\times 10^{-4}$
$5.381\times 10^{-4}$ $6.669\times 10^{-4}$ $92$ $3.799\times 10^{-6}$
$2.149\times 10^{-5}$ $7.418\times 10^{-5}$ $1.212\times 10^{-4}$ $1.525\times
10^{-4}$ $1.570\times 10^{-4}$ $1.456\times 10^{-4}$ $1.061\times 10^{-4}$
$93$ $3.383\times 10^{-6}$ $1.659\times 10^{-5}$ $4.472\times 10^{-5}$
$5.729\times 10^{-5}$ $5.323\times 10^{-5}$ $4.434\times 10^{-5}$ $3.023\times
10^{-5}$ $1.330\times 10^{-5}$ $94$ $2.670\times 10^{-6}$ $1.138\times
10^{-5}$ $2.408\times 10^{-5}$ $2.435\times 10^{-5}$ $1.686\times 10^{-5}$
$1.145\times 10^{-5}$ $5.806\times 10^{-6}$ $1.571\times 10^{-6}$ $95$
$1.978\times 10^{-6}$ $7.340\times 10^{-6}$ $1.226\times 10^{-5}$ $9.851\times
10^{-6}$ $5.124\times 10^{-6}$ $2.855\times 10^{-6}$ $1.084\times 10^{-6}$
$1.825\times 10^{-7}$ $96$ $1.411\times 10^{-6}$ $4.571\times 10^{-6}$
$6.056\times 10^{-6}$ $3.883\times 10^{-6}$ $1.527\times 10^{-6}$ $7.009\times
10^{-7}$ $2.006\times 10^{-7}$ $2.114\times 10^{-8}$ $97$ $9.818\times
10^{-7}$ $2.784\times 10^{-6}$ $2.938\times 10^{-6}$ $1.510\times 10^{-6}$
$4.515\times 10^{-7}$ $1.712\times 10^{-7}$ $3.704\times 10^{-8}$ $2.453\times
10^{-9}$ $98$ $6.709\times 10^{-7}$ $1.667\times 10^{-6}$ $1.405\times
10^{-6}$ $5.810\times 10^{-7}$ $1.324\times 10^{-7}$ $4.150\times 10^{-8}$
$6.798\times 10^{-9}$ $2.830\times 10^{-10}$

By using the present cosmogenic neutrino flux upper limit [30, 31, 32],
$\Phi(E)$, we evaluate the expected number of $\tau$ events in the TAUWER
array according to the following formula

$N_{\tau}=\int_{E_{min}}^{E_{max}}\Phi(E)P_{escp}(E)(\varepsilon_{acc}(E,\Theta,\Phi,L^{*})\times\beta_{\tau})dE\Delta
t\Delta A\Delta\Omega$

(3)

where $\beta_{\tau}$ is $\tau$ branching ratios for $\pi^{-}$,
$\pi^{-}\pi^{0}$, $\pi^{-}\pi^{+}\pi^{-}$, $\pi^{-}\pi^{0}\pi^{0}$ and
$\pi^{-}\pi^{+}\pi^{-}\pi^{0}$ decay modes, $\varepsilon_{acc}$ is the
efficiency of detector array acceptance. $\Delta t$ is the time interval of
one year, $\Delta A=0.6$ $km^{2}$ area and $\Delta\Omega$ is the solid angle
that watched the showers emerging below the horizon within $8\degree$ in
$\theta$ and the showers emerging from the chain of mountain in front of the
array of $2\degree$ in $\theta$ . To evaluate the number of detected events,
we take into account two scenarios: Earth Skimming and the screen strategies
produced by a mountain chain in front of the array.

The expected number of events by the TAUWER array was estimated using the
latest neutrino flux upper limit reported by IceCube experiment [30, 31, 32].
The results of different scenarios, only skimming strategy and screen strategy
are summarized in Table 2. The expected number of events as function of two
different neutrino flux upper limits are $0.18$ and $0.43$ events in $km^{2}$
per year for skimming and skimming-screen strategy, respectively. If the
integrated time is 3 years [32] and the surface of array is 10 $km{{}^{2}}$
the number of expected events for TAUWER array are about $13$ events. The
presence of ‘screen’ increases the detected number of events about a factor
three. That is explained by higher probability of interaction of the
$\nu_{\tau}$ that travels a longer path when it is skimmed by the Earth crust.

Table 2: The expected number of events in $km^{2}$ per year for TAUWER array
using analytic method and DANTON simulation code corresponding to the
different upper limit fluxes. With DANTON simulation we computed only skimming
events.

Flux Scenario Energy range $N_{event}$ $N_{event}$ ($GeVcm^{-2}s^{-1}sr^{-1}$)
(PeV) CTEQ6.6M DantonSim Skimming $[10-1000]$ 0.0437 0.0885 $2\times 10^{-8}$
(all flavor) [30] Screen $[10-1000]$ 0.1273 Skimming + Screen $[10-1000]$
0.1710 Skimming $[10-1000]$ 0.1116 0.2257 $5.1\times 10^{-8}$ ($\nu_{\tau}$
flavor) [32] Screen $[10-1000]$ 0.3245 Skimming + Screen $[10-1000]$ 0.4361

The $P_{escp}$ computed using the analytic method has been compared to that
evaluated using a specific code DANTON that simulates the tau-neutrino
interaction in the Earth crust and the $\tau$ decay giving the $\tau$ decay
point in the air [33]. Using this code we obtained the total $P_{escp}$ of
$\tau$ for energy greater than 200 PeV in average a factor 2.5 higher than the
$P_{escp}$ evaluated with the analytic method. This discrepancy is due mainly
to the different PDF used in the neutrino cross section parametrization. To
evaluate the expected $\tau$ lepton we can use the geometrical acceptance
computed in our study because the tau decay point coordinates in DANTON have
an average distance from the array of 8 km and a momentum distribution similar
to which was used in the analytic method. By using the upper limit fluxes [30,
32], we obtain 0.21 and 0.55 events per km2 per year with the DANTON
simulation code.

###  3.2 Arrival $\tau$ shower Direction

Figure 4:  Zenith (top) and Azimuthal (bottom) angles evaluated from active
station for different showers where decay length from emerging point is 7 km
and primary energy ranges from $10$ PeV to $1000$ PeV. Active stations are
selected by requiring at least two hits on the tower.

The direction of the $\tau$ air shower is calculated with a common method used
in literature [34] based on the minimization of the direction vectors of each
track by using Minuit [35]. The momentum of each particle are used to obtain
the direction vectors $(n_{x}$, $n_{y}$, $n_{z})$ of the shower in Eq. 4.

$\displaystyle\left.\begin{aligned}
\chi^{2}=\sum_{i}^{N}w_{i}\Big{(}&n_{x}(x_{i}-x_{core})+n_{y}(y_{i}-y_{core})\\\
&+n_{z}(z_{i}-z_{core})-c(t_{i}-t_{core})\Big{)}^{2}\end{aligned}\right.$ (4)

where $t_{i}$ is the arrival time of $e^{\pm}$ , $x$, $y$ and $z$ its
coordinates on the array plane, $x_{core},\,y_{core},\,z_{core}$ the
coordinates of the core of the shower evaluated with maximum weighted density
of the particles on the inclined plane. $t_{core}$ is estimated from the
closest station taken into account of the distance from the core. The active
stations are selected by requiring at least two particles on the tower and
used in the evaluation of the shower direction for different showers where the
primary in an energy range of $10$ PeV - $1000$ PeV and the decay point 7 km
distant from the array plane.

Figure 4 shows both zenith and azimuthal angles evaluated for sample of $27$
showers in the energy range between 10 and 1000 PeV. The arrival directions
$\theta$ and $\phi$ are estimated with an error of $0.16\degree$ and
$0.045\degree$, respectively.

###  3.3 $\tau$ shower energy reconstruction

Figure 5: a) Number of hits per station versus distance from the shower core
for different $\tau$ energies, b) the energy of $\tau$ showers as function of
the shower radius.

The center of the shower is estimated from the most triggered stations and
hits per station as function of distance from the shower core is given in
Figure 5. This figure shows the number of hits per station as a function of
distance for different energy intervals between 10 and 1000 PeV where at least
one hit on a single tile is required. The core of the shower is evaluated from
the stations having at least 5 hits. While the radius of the shower is
obtained by requiring 4 subsequent stations having hits ranging from 1 to 4
and register the first one as the radius (edge) of the shower. The Figure 5
shows the energy versus its corresponding radius. It is clearly seen that
there is an exponential behaviour between energy and radius in the range of
the interested energy region.

###  3.4 Trigger Decision

The average timing spectrum for $e^{\pm}$ between 10 PeV and 1000 PeV for all
hadronic channels of $\tau$ shower that hit any scintillator in the 640
stations has been studied. The timing spread of the stations is computed by
using the arrival time of all hits in the station and identifying the earliest
one which is selected as time reference to find a delay time in a $(4\times
4)$ subarray around the center of the shower. The maximum time difference
between the first and last hit produced by the $\tau$ shower on the array
plane must be in the time interval of $0.3\,\mu s$. From these simulation
studies the signature of $\tau$ shower is given by a cluster of station with
highest density of $e^{\pm}$ and the maximum number of hits in a single
station for $\tau$ shower between 10 PeV and 1000 PeV ranges from $450-3500$,
in the time interval of $10-20$ ns. Compared this requirement to the vertical
EAS events collected in a test made at Karlsruhe Institute of Technology (KIT)
[19], the atmospheric background can be easily suppressed if we require $250$
hits, defined as trigger T1. The second level trigger T2 can be defined as
requiring at least $450$ hits on each station of a sub-array $(2\times 2)$ in
a time window of $88$ ns. This criteria (/trigger) leads to a selection of a
shower with an energy of $10^{8}\,GeV$ or greater.

##  4 Conclusions

We have simulated $\tau$ decay modes ($\pi^{-}$, $\pi^{-}\pi^{0}$,
$\pi^{-}\pi^{+}\pi^{-}$, $\pi^{-}\pi^{0}\pi^{0}$,
$\pi^{-}\pi^{+}\pi^{-}\pi^{0}$) that covers 64 % of the $\tau$ branching ratio
in the energy range of $10-1000$ PeV to determine the performance of the
TAUWER array. It shows that the energy of the showers between $100$ and $1000$
PeV can be reconstructed in terms of the radius as discussed in Section 3.3.
The arrival direction is also determined with an accuracy of $0.16^{o}$ in the
same energy range. The atmospheric background can be easily rejected by
requiring T1 trigger while a second level trigger, T2 selects showers with
$10^{8}$ GeV or greater energies. We have estimated the $\nu_{\tau}$
conversion to $\tau$ and $\tau$ propagation through the rock including both
skimming and screen strategies with the TAUWER detector array, which is
ranging between $5$-$13$ tau lepton events in the energy range of $10-1000$
PeV for $3$ years of integrated time period and an array surface of 10
$km^{2}$ using the analytical method. We obtained $6$-$16$ tau lepton skimmed
events by using DANTON simulation for the same time period, array surface and
in same energy range. The number of events skimmed by the Earth crust using
this simulation code is higher of a factor two due to the neutrino cross
section. The effect of the screen in front of the array increase the expected
number of events of a factor three.

##  5 Acknowledgement

We acknowledge the financial support by the Scientific and Technological
Research Council of Turkey (TUBITAK) (grant number: 114F138).

## References

  * [1] F. Halzen and E. Zas, “Giant horizontal air showers. Implications for AGN neutrino fluxes,” Physics Letters B, vol. 289, no. 1, pp. 184–188, 1992\.
  * [2] J. Linsley and L. Scarsi, “Cosmic-Ray Composition at ${10}^{17}$ \- ${10}^{18}$ eV,” Phys. Rev. Lett., vol. 9, pp. 123–125, Aug 1962.
  * [3] F. Kakimoto, S. Ogio, Y. Shirasaki, N. Gironda, T. Kaneko, H. Yoshii, T. Baba, K. Nishi, N. Tajima, Y. Yamada, S. Shimoda, H. Nakatani, E. Gotoh, Y. Toyoda, K. Murakami, Y. Matsubara, Y. Mizumoto, W. Velasquez, and A. Velarde, “A small air-shower array at Mount Chacaltaya,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 373, no. 2, pp. 282–289, 1996.
  * [4] S. Barwick et al., “A first search for cosmogenic neutrinos with the ARIANNA Hexagonal Radio Array,” Astroparticle Physics, vol. 70, pp. 12–26, 2015.
  * [5] M. Ackermann, J. Adams, J. Ahrens, K. Andeen, J. Auffenberg, X. Bai, B. Baret, S. W. Barwick, R. Bay, K. Beattie, and et al., “Search for Ultra–High‐Energy Neutrinos with AMANDA‐II,” The Astrophysical Journal, vol. 675, pp. 1014–1024, Mar 2008.
  * [6] P. Sokolsky, “Final Results from the High solution Fly’s Eye (HiRes) Experiment,” Nuclear Physics B - Proceedings Supplements, vol. 212-213, pp. 74–78, 2011.  Proceedings of the Cosmic Ray International Seminars (CRIS 2010).
  * [7] R. Abbasi et al., “Search for ultrahigh-energy tau neutrinos with IceCube,” Phys. Rev. D, vol. 86, p. 022005, Jul 2012.
  * [8] P. Abreu et al., “A search for ultra-high energy neutrinos in highly inclined events at the Pierre Auger Observatory,” Phys. Rev. D 84, 122005 (2011); Erratum: Phys. Rev. D 85, 029902(E) (2012), Feb. 2012.
  * [9] IceCube Collaboration, “Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector,” Science, vol. 342, pp. 1242856–1242856, Nov 2013.
  * [10] IceCube Collaboration, Fermi-LAT, MAGIC, AGILE, ASAS-SN, HAWC, H.E.S.S., INTEGRAL, Kanata, Kiso, Kapteyn, L. Telescope, Subaru, Swift/NuSTAR, and VERITAS teams, “Multimessenger observations of a flaring blazar coincident with high-energy neutrino IceCube-170922A,” Science, vol. 361, Jul 2018.
  * [11] E. Zas, “Neutrino Detection with Inclined Air Showers,” New Journal of Physics, vol. 7, Apr. 2005.
  * [12] D. Fargion, A. Aiello, and R. Conversano, “Horizontal tau air showers from mountains in deep valley: Traces of UHECR neutrino tau,” in Proceedings, 26th International Cosmic Ray Conference (ICRC), August 17-25, 1999, Salt Lake City: Invited, Rapporteur, and Highlight Papers, p. 396, 1999\.  [2,396(1999)].
  * [13] Y. Asaoka and M. Sasaki, “Cherenkov $\tau$ shower earth-skimming method for PeV-EeV $\nu_{\tau}$ observation with Ashra,” Astroparticle Physics, vol. 41, pp. 7–16, 2013.
  * [14] D. Gra, M. Manganaro, E. Bernardini, M. Doro, M. Will, S. Lombardi, J. Rico, and D. Sobczynska, “Search for tau neutrinos at PeV energies and beyond with the MAGIC telescopes,” PoS, vol. EPS-HEP2017, p. 017, 2017.
  * [15] Z. Cao, M. Huang, P. Sokolsky, and Y. Hu, “Ultra high energy $\nu_{\tau}$ detection with a cosmic ray tau neutrino telescope using fluorescence/Cerenkov light technique,” Journal of Physics G: Nuclear and Particle Physics, vol. 31, no. 7, p. 571, 2005.
  * [16] M. Meier, J. Soedingrekso, and T. IceCube.Collaboration, “Search for astrophysical tau neutrinos with an improved double pulse method,” Proceedings of Science, vol. 358, p. 960, July 2019.
  * [17] D. R. Williams and T. IceCube.Collaboration, “Results from icecube,” Proceedings of Science, vol. 358, p. 016, July 2019.
  * [18] J. Nam and T. Liu, “Feasibility of antenna array experiment for Earth skimming tau-neutrino detection in Antarctica ,” Proceedings of Science, vol. 301, Aug. 2018.
  * [19] M. Iori, E. Arslan, H. Denizli, M. Kaya, A. Yilmaz, and J. Russ, “Electron–muon identification by atmospheric shower in a new concept of an EAS detector,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 692, pp. 285–287, 2012.  3rd Roma International Conference on Astroparticle Physics.
  * [20] M. Iori, H. Denizli, A. Yilmaz, F. Ferrarotto, and J. Russ, “Electron–muon identification by atmospheric shower and electron beam in a new eas detector concept,” The Astrophysical Journal, vol. 801, p. 140, Mar 2015.
  * [21] S. Jadach, Z. Was, R. Decker, and J. Kühn, “The tau decay library TAUOLA: Version 2.4,” Computer Physics Communications, vol. 76, no. 3, pp. 361 – 380, 1993.
  * [22] D. Heck, G. Schatz, T. Thouw, J. Knapp, and J. N. Capdevielle, “CORSIKA: A Monte Carlo code to simulate extensive air showers,” 1998.
  * [23] S. Ostapchenko”, “Nonlinear screening effects in high energy hadronic interactions,” Physical Review D, vol. 74, no. 1, 2006.
  * [24] H. Fesefeldt, “The simulation of hadronic showers: physics and applications,” tech. rep., 1985.
  * [25] Y. S. Jeong, M. V. Luu, M. H. Reno, and I. Sarcevic, “Tau energy loss and ultrahigh energy skimming tau neutrinos,” Phys. Rev. D, vol. 96, p. 043003, Aug 2017.
  * [26] R. Gandhi, C. Quigg, M. H. Reno, and I. Sarcevic, “Neutrino interactions at ultrahigh energies,” Phys. Rev. D, vol. 58, p. 093009, Sep 1998.
  * [27] Y. S. Jeong and M. H. Reno, “Quark mass effects in high energy neutrino-nucleon scattering,” Phys. Rev. D, vol. 81, p. 114012, Jun 2010\.
  * [28] M. H. Reno, T. M. Venters, J. F. Krizmanic, L. A. Anchordoqui, C. Guepin, and A. V. Olinto, “A new calculation of Earth-skimming very- and ultra-high energy tau neutrinos,” PoS, vol. ICRC2019, p. 989, 2019.
  * [29] S. I. Dutta”, “Tau neutrino propagation and tau energy loss,” Physical Review D, vol. 72, no. 1, 2005.
  * [30] M. G. Aartsen, M. Ackermann, J. Adams, et al., “Differential limit on the extremely-high-energy cosmic neutrino flux in the presence of astrophysical background from nine Years of IceCube data,” Phys. Rev. D, vol. 98, p. 062003, Sep 2018.
  * [31] M. G. Aartsen, M. Ackermann, et al., “Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data,” Phys. Rev. Lett., vol. 113, p. 101101, Sep 2014.
  * [32] M. G. Aartsen et al., “Search for astrophysical tau neutrinos in three Years of IceCube data,” Physical Review D, vol. 93, no. 2, 2016.
  * [33] V. Niess and O. Martineau-Huynh, “DANTON: a Monte-Carlo sampler of $\tau$ from $\nu_{\tau}$ interacting with the Earth,” 2018.
  * [34] H. H. Kh, “A new model-independent approach for finding the arrival direction of an extensive air shower,” Journal of Cosmology and Astroparticle Physics, Nov. 2016.
  * [35] F. James and M. Roos, “Minuit - a system for function minimization and analysis of the parameter errors and correlations,” Computer Physics Communications, vol. 10, no. 6, pp. 343–367, 1975.