# Large Deviations of Telecom Processes

M.A. Lifshits  S.E. Nikitin

(July 2021)

###### Abstract

We study large deviation properties of Telecom processes appearing as limits
in a critical regime of infinite source Poisson models.

AMS Subject Classification: Primary: 60F10. Secondary: 60F05, 60G57, 60E07.

Keywords: large deviation probabilities, Telecom process, Poisson random
measure, teletraffic, workload.

##  1 Introduction: Telecom processes

###  1.1 A service system

Telecom processes originate from a remarkable work by I. Kaj and M.S. Taqqu
[9] who handled the limit behavior of ”teletraffic systems” by using the
language of integral representations as a unifying technique. Their article
brightly represents a wave of interest to the subject, see e.g. [8, 10, 12,
13, 14], and the surveys with further references [5, 6, 7], to mention just a
few. Simplicity of the dependence mechanism used in the model enables to get a
clear understanding both of long range dependence in one case, and independent
increments, in other cases.

The work of the system represents a collection of _service processes_ or
_sessions_ , using telecommunication terminology. Every process starts at some
time $s$, lasts $u$ units of time, and occupies $r$ _resource_ units (synonyms
for resource are _reward, transmission rate_ etc). The amount of occupied
resources $r$ remains constant during every service process.

The formal model of the service system is based on Poisson random measures and
looks as follows. Let
$\mathcal{R}:=\\{(s,u,r)\\}={\mathbb{R}}\times{\mathbb{R}}_{+}\times{\mathbb{R}}_{+}$.
Every point $(s,u,r)$ corresponds to a possible service process with starting
time $s$, duration $u$, and required resources $r$.

The system is characterized by the following parameters:

  * •

$\lambda>0$ – _arrival intensity_ of service processes;

  * •

$F_{U}(du)$ – the distribution of service duration;

  * •

$F_{R}(dr)$ – the distribution of amount of required resources.

One may assume ${\mathbb{P}}(R>0)=P(U>0)=1$ without loss of generality.

Define on $\mathcal{R}$ an intensity measure

$\mu(ds,du,dr)=\lambda ds\,F_{U}(du)\,F_{R}(dr).$

Let $N$ be a Poisson random measure with intensity $\mu$. One may consider the
samples of $N$ (sets of triplets $(s,u,r)$, each triplet corresponding to a
service process) as variants (sample paths) of the work for the system.

The instant workload on the system at time $t$ writes as

$W^{\circ}(t)=\int_{\mathcal{R}}r{\mathbf{1}}_{\\{s\leq t\leq s+u\\}}dN.$

This is essentially the sum of occupied resources over the processes active at
time $t$. The integral workload over the interval $[0,t]$ is

$\displaystyle W^{*}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}W^{\circ}(\tau)d\tau=\int_{\mathcal{R}}r\int_{0}^{t}{\mathbf{1}}_{\\{s\leq\tau\leq
s+u\\}}d\tau dN$ $\displaystyle=$
$\displaystyle\int_{\mathcal{R}}r\cdot\Big{|}[s,s+u]\cap[0,t]\Big{|}dN:=\int_{\mathcal{R}}r\ell_{t}(s,u)dN.$

Here $|\cdot|$ stands for the length of an interval, and the kernel

$\ell_{t}(s,u):=\Big{|}[s,s+u]\cap[0,t]\Big{|}$ (1)

will be often used in the sequel.

Notice that $W^{\circ}(\cdot)$ is a stationary process and its integral
$W^{*}(\cdot)$ is a process with stationary increments.

We suppose that either the variables $R$ and $U$ have finite variance, or
their distributions have regular tails. More precisely, either

${\mathbb{P}}(U>u)\sim\frac{{c_{\scriptscriptstyle U}}}{u^{\gamma}}\ ,\qquad
u\to\infty,\qquad 1<\gamma<2,\ {c_{\scriptscriptstyle U}}>0,$

or ${\mathbb{E}}U^{2}\,<\infty$. In the latter case we formally set
$\gamma:=2$.

Analogously, we assume either

${\mathbb{P}}(R>r)\sim\frac{{c_{\scriptscriptstyle R}}}{r^{\delta}}\ ,\qquad
r\to\infty,\qquad 1<\delta<2,\ {c_{\scriptscriptstyle R}}>0,$

or ${\mathbb{E}}R^{2}\,<\infty$. In the latter case we formally set
$\delta:=2$.

The behavior of the service system crucially depends of the parameters
$\gamma,\delta\in(1,2]$.

###  1.2 Limit theorems for the workload

####  1.2.1 Centered and scaled workload process

The main object of theoretical interest is the behavior of the integral
workload as a process (function of time) observed on long time intervals.

In order to obtain a meaningful limit, one must

scale (contract) the time so that it would run through the standard time
interval, center the workload process, and divide it by an appropriate scalar
factor.

We choose $[0,1]$ as a standard time interval.

Centering and scaling by appropriate factor $b$ lead to a _normalized integral
workload process_

$Z_{a}(t):=\frac{W^{*}(at)-{\mathbb{E}}R\,\cdot{\mathbb{E}}\,U\,\cdot a\lambda
t}{b}\ ,\qquad t\in[0,1],\ b=b(a,\lambda).$

####  1.2.2 A limit theorem leading to a Telecom process

It is remarkable that a simple tuning of three parameters
$\lambda,\gamma,\delta$ may lead to different limiting processes for $Z_{a}$,
namely, one can obtain

  * •

a Wiener process;

  * •

a fractional Brownian motion with index $H\in(1/2,1)$;

  * •

a centered Lévy stable process with positive spectrum;

  * •

a stable Telecom process;

  * •

a Poisson Telecom process.

While the first three processes present a core of the classical theory of
stochastic processes, the Telecom processes are almost not studied. In this
article we focus on some key properties of the Poisson Telecom process.

For the full panorama of related limit theorems we refer to [11, Chapter 3]
and recall here only one result concerning the Poisson Telecom process (cf.
[11, Theorem 13.16]) related to the case of _critical intensity_

$\frac{\lambda}{a^{\gamma-1}}\to L,\qquad\qquad 0<L<\infty.$ (2)

###### Theorem 1

Assume that $1<\gamma<\delta\leq 2$, $a\to\infty$, and that critical intensity
condition $\eqref{critint}$ holds. Let $Q:=L\,{c_{\scriptscriptstyle
U}}\,\gamma$. Then with scaling $b:=a$ the finite dimensional distributions of
the process $Z_{a}$ converge to those of the Poisson Telecom process
$Y_{Q,\gamma}$ admitting an integral representation

$Y_{Q,\gamma}(t)=\int_{\mathcal{R}}r\,\ell_{t}(s,u)\bar{N}_{Q,\gamma}(ds,du,dr).$

Here $\ell_{t}(s,u)$ is the kernel defined in $(\ref{ell})$ and
$\bar{N}_{Q,\gamma}$ is a centered Poisson random measure of intensity
$Q\,\mu_{\gamma}$ where

$\mu_{\gamma}(ds,du,dr):=\frac{ds\,du}{u^{\gamma+1}}\ F_{R}(dr).$

For studies on Poisson Telecom process we refer to [2, 4]. It is well known
that the process $(Y_{Q,\gamma}(t))_{t\geq 0}$, is correctly defined if
${\mathbb{E}\,}(R^{\gamma})<\infty$. In accordance with its role in the limit
theorem, it has stationary increments. It is, however, not self-similar like
other limiting processes in the same model, such as Wiener process, fractional
Brownian motion, or strictly stable Lévy process.

##  2 Main results

###  2.1 A limit theorem for Telecom process

At large time scales the Poisson Telecom process essentially behaves as a
$\gamma$-stable Lévy process. This fact is basically known but we present it
here for completeness of exposition. The analogy with a stable law will also
guide us (to some extent and within a certain range) in the subsequent studies
of large deviation probabilities.

###### Proposition 2

We have a weak convergence

$\left({\mathbb{E}\,}(R^{\gamma})\,t\right)^{-1/\gamma}\,Y_{Q,\gamma}(t)\
\Rightarrow\SS_{Q,\gamma},\qquad\textrm{ as }t\to\infty,$ (3)

where $\SS_{Q,\gamma}$ is a centered strictly $\gamma$-stable random variable
with positive spectrum, i.e.

${\mathbb{E}\,}\exp\\{it\SS_{Q,\gamma}\\}=\exp\left\\{Q\int_{0}^{\infty}\frac{e^{itu}-1-itu}{u^{\gamma+1}}\,du\right\\}$

###  2.2 Large deviations

According to the limit theorem (2), large deviation probability is

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho),\quad\textrm{as }\
\varrho=\varrho(t)\gg t^{1/\gamma}.$

Its behavior may be different in different zones of $\varrho$ and may depend
on the distribution of $R$. We consider the most important cases in the
following subsections.

####  2.2.1 Moderate large deviations

###### Theorem 3

Let $\varrho=\varrho(t)$ be such that $t^{1/\gamma}\ll\varrho\ll t$. Then

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)=D\,t\,\varrho^{-\gamma}\
(1+o(1)),\qquad\textrm{as }\ t\to\infty,$ (4)

where $D:=\tfrac{Q\ {\mathbb{E}\,}(R^{\gamma})}{\gamma}$.

This result should be compared with the limit theorem (3) because (4) yields

$\displaystyle{\mathbb{P}}\left(({\mathbb{E}\,}(R^{\gamma})t)^{-1/\gamma}Y_{Q,\gamma}(t)\geq\rho\right)={\mathbb{P}}\left(Y_{Q,\gamma}(t)\geq({\mathbb{E}\,}(R^{\gamma})t)^{1/\gamma}\rho\right)$
$\displaystyle\sim$ $\displaystyle
D\,t\,({\mathbb{E}\,}(R^{\gamma})t)^{-1}\rho^{-\gamma}=\frac{Q}{\gamma}\
\rho^{-\gamma}\sim{\mathbb{P}}(\SS_{Q,\gamma}\geq\rho),$

whenever $1\ll\rho\ll t^{-(\gamma-1)/\gamma}$. In other words, the moderate
large deviation probabilities are equivalent to those of the limiting
distribution.

Using the terminology of the background service system, moderate deviation is
attained by a unique heavy service process. We will stress this fact later in
the proof.

####  2.2.2 Intermediate large deviations

The following result describes the situation on the upper boundary of moderate
deviations’ zone.

###### Theorem 4

Let $\kappa>0$ be such that ${\mathbb{P}}(R\geq\kappa)>0$ and

${\mathbb{P}}(R=\kappa)=0.$ (5)

Let $\varrho=\varrho(t)=\kappa t$. Then

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)=QD^{(1)}_{I}(\kappa)\,t^{-(\gamma-1)}\,(1+o(1)),\qquad\textrm{as
}\ t\to\infty,$

where

$D_{I}^{(1)}(\kappa):=\left(\frac{\kappa^{-\gamma}}{\gamma}\,{\mathbb{E}\,}(R^{\gamma}{\mathbf{1}}_{\\{R\geq\kappa\\}})+\frac{(2-\gamma)\kappa^{1-\gamma}}{(\gamma-1)\gamma}\,{\mathbb{E}\,}(R^{\gamma-1}{\mathbf{1}}_{\\{R\geq\kappa\\}})\right).$

###### Remark 5

There is a continuity between the moderate and intermediate zones in what
concerns the degree of $t$ but the constant in the intermediate case is
different. Indeed, by plugging formally $\rho:=\kappa t$ into (4) one obtains
the asymptotics
$\tfrac{Q\,{\mathbb{E}\,}(R^{\gamma})}{\gamma}\,\kappa^{-\gamma}\,t^{-(\gamma-1)}$
which corresponds to the fist term in the definition of $D_{I}^{(1)}(\kappa)$.
When $\kappa$ goes to zero, the second term in that definition is smaller than
the first one because

$\kappa\,{\mathbb{E}\,}(R^{\gamma-1}{\mathbf{1}}_{\\{R\geq\kappa\\}})\leq{\mathbb{E}\,}(R^{\gamma}{\mathbf{1}}_{\\{R\geq\kappa\\}})\searrow
0,\qquad\textrm{as }\kappa\to 0.$

###### Remark 6

If (5) does not hold, the decay order of large deviations will be the same but
the expression for the corresponding constant becomes more involved and less
explicit.

The attentive reader will notice that Theorem 4 does not work for large
$\kappa$ if the distribution of $R$ is compactly supported. Indeed, in this
case the large deviation asymptotics will be different, as the next result
shows. In terms of the service system, it handles the case when the large
deviation can be attained by accumulation of $n$ heavy service processes but
cannot be attained by $(n-1)$ ones.

###### Theorem 7

Let $\kappa>0$. Let $n$ be the positive integer such that
${\mathbb{P}}(R\geq\tfrac{\kappa}{n})>0$ but

${\mathbb{P}}(R\geq\tfrac{\kappa}{n-\zeta})=0\qquad\textrm{for some
}\zeta\in(0,1).$ (6)

Assume that

${\mathbb{P}}(R_{1}+\cdots+R_{n}=\kappa)=0,$ (7)

where $R_{1},\dots,R_{n}$ are independent copies of $R$.

Let $\varrho=\varrho(t):=\kappa t$. Then

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)=Q^{n}D^{(n)}_{I}(\kappa)\,t^{-(\gamma-1)n}\,(1+o(1)),\qquad\textrm{as
}\ t\to\infty,$

where $D^{(n)}_{I}(\kappa)$ is some finite positive constant depending on
$n,\kappa$ and on the law of $R$.

###### Remark 8

The explicit form of $D^{(n)}_{I}(\kappa)$ is given in equation (31) below.

###### Remark 9

Theorem 7 does not cover a critical case $\zeta=1$, where we have
${\mathbb{P}}(R\geq\tfrac{\kappa}{n-1})=0$ but
${\mathbb{P}}(R\geq\tfrac{\kappa}{n-1}-\varepsilon)>0$ for all
$\varepsilon>0$. In this case, the assertion of the theorem may not hold
because the large deviation probability behavior depends of that of the upper
tail
${\mathbb{P}}(R\in[\tfrac{\kappa}{n-1}-\varepsilon,\tfrac{\kappa}{n-1}))$, as
$\varepsilon\to 0$.

####  2.2.3 Ultralarge deviations

###### Theorem 10

Let $\varrho=\varrho(t)\gg t$. Assume that the tail probability function
${\bar{F}_{R}}(y):={\mathbb{P}}(R\geq y)$ is regularly varying of negative
order $-m$ where $m>\gamma$. Then

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)=Q\,D\,t^{-(\gamma-1)}{\bar{F}_{R}}(\varrho/t)\,(1+o(1)),\qquad\textrm{as
}\ t\to\infty,$

where

$D:=\frac{m(m-1)}{\gamma(\gamma-1)(m-\gamma+1)(m-\gamma)}.$

As in Theorem 3, the workload’s large deviation is attained by a unique long
and heavy service process.

Theorem 10 deals with the distributions of $R$ having essentially polynomial
tails. The corresponding distributions with light tails lead to completely
different results such as Poisson large deviations. This direction requires
supplementary research to be presented elsewhere.

###  2.3 Concluding remark

A challenging case when the workloads’ ultralarge deviation is formed via the
interaction of infinitely many service processes remains beyond the scope of
this article. Here, a large deviation rate function related to the
distribution of $R$ must play a major role and the results in the spirit of
classical large deviation theory [3] are expected. This might be a subject of
a subsequent work.

##  3 Proofs

###  3.1 Preliminaries

Let us introduce two auxiliary intensity measures. The first one is the
”distribution” of the kernel $\ell_{t}$, namely

$\mu^{(\ell)}_{t}(A):=\int_{{\mathbb{R}}}\int_{{\mathbb{R}}_{+}}{\mathbf{1}}_{\\{\ell_{t}(s,u)\in
A\\}}\frac{du}{u^{\gamma+1}}\,ds,\qquad A\in\mathcal{B}([0,t]).$

The second is the ”distribution” of the product $r\ell_{t}(s,u)$,

$\displaystyle\mu^{(\ell,r)}_{t}(A)$ $\displaystyle:=$
$\displaystyle\int_{{\mathbb{R}}_{+}}\int_{{\mathbb{R}}_{+}}{\mathbf{1}}_{\\{r\ell\in
A\\}}\mu^{(\ell)}_{t}(d\ell)F_{R}(dr)$ $\displaystyle=$
$\displaystyle\mu\\{(s,u,r):r\ell_{t}(s,u)\in A\\},\qquad
A\in\mathcal{B}({\mathbb{R}}_{+}).$

A simple variable change in the definition of $Y_{Q,\gamma}(t)$ yields

$Y_{Q,\gamma}(t)=\int_{{\mathbb{R}}_{+}}v{\widetilde{N}}_{Q,\gamma}(dv)$ (8)

where ${\widetilde{N}}_{Q,\gamma}$ is a centered Poisson measure with
intensity $Q\mu^{(\ell,r)}_{t}$. Therefore, the properties of
$\mu^{(\ell,r)}_{t}$ determine those of $Y_{Q,\gamma}(t)$.

As a first step, we give an explicit formula for the intermediate measure
$\mu^{(\ell)}_{t}$. First, by definition we have
$\mu^{(\ell)}_{t}(t,\infty)=0$. Next, let us fix an $\ell_{0}\in(0,t]$ and
find $\mu^{(\ell)}_{t}[\ell_{0},t]$. In fact, $\ell_{t}(s,u)\geq\ell_{0}$ iff
$u\geq\ell_{0}$ and $s\in[\ell_{0}-u,t-\ell_{0}]$. Therefore,

$\mu^{(\ell)}_{t}[\ell_{0},t]=\int_{\ell_{0}}^{\infty}(t-2\ell_{0}+u)\frac{du}{u^{\gamma+1}}=\frac{t\,\ell_{0}^{-\gamma}}{\gamma}+\frac{2-\gamma}{(\gamma-1)\gamma}\,\ell_{0}^{1-\gamma}.$
(9)

It follows that the measure $\mu^{(\ell)}_{t}$ has a weight
$\frac{t^{-(\gamma-1)}}{(\gamma-1)\gamma}$ at the right boundary point $t$ and
a density

$\frac{d\mu^{(\ell)}_{t}}{d\ell}(\ell)=t\,\ell^{-1-\gamma}+\frac{2-\gamma}{\gamma}\,\ell^{-\gamma},\qquad
0<\ell<t.$

For each $\ell_{0}>0$, formula (9) also yields a bound

$\mu^{(\ell)}_{t}[\ell_{0},\infty)=\mu^{(\ell)}_{t}[\ell_{0},t]\leq\frac{t\ell_{0}^{-\gamma}}{\gamma}\left(1+\frac{2-\gamma}{\gamma-1}\right)=\frac{t\,\ell_{0}^{-\gamma}}{\gamma(\gamma-1)}.$
(10)

Finally, consider the asymptotic behavior of

$\mu^{(\ell,r)}_{t}[\varrho,\infty)=\int_{{\mathbb{R}}_{+}}\mu^{(\ell)}_{t}\left[\frac{\varrho}{r},t\right]F_{R}(dr).$
(11)

Assume that $\varrho\to\infty$ but $\varrho/t\to 0$. Then it follows from (9)
that for every fixed $r$

$\mu^{(\ell)}_{t}\left[\frac{\varrho}{r},t\right]=\frac{t\,\varrho^{-\gamma}r^{\gamma}}{\gamma}(1+o(1)).$
(12)

By using (10), we also have an integrable majorant w.r.t. the law $F_{R}$:

$\mu^{(\ell)}_{t}\left[\frac{\varrho}{r},t\right]\leq\frac{t\,\varrho^{-\gamma}r^{\gamma}}{\gamma(\gamma-1)}.$

By integrating this estimate in (11) we obtain

$\mu^{(\ell,r)}_{t}[\varrho,\infty)\leq\frac{{\mathbb{E}\,}(R^{\gamma})}{\gamma(\gamma-1)}\,t\,\varrho^{-\gamma}.$
(13)

Furthermore, by Lebesgue’s majorated convergence theorem (11) and (12) yield

$\mu^{(\ell,r)}_{t}[\varrho,\infty)=\varrho^{-\gamma}\,t\int_{{\mathbb{R}}_{+}}\frac{r^{\gamma}}{\gamma}F_{R}(dr)\,(1+o(1))=\frac{{\mathbb{E}\,}(R^{\gamma})}{\gamma}\,t\,\varrho^{-\gamma}\,(1+o(1)).$
(14)

###  3.2 Proof of Proposition 2

Consider the integral representation (8). According to a general criterion of
the weak convergence of Poisson integrals to a stable law [11, Corollary 8.5],
it is enough to check that for each fixed $\rho>0$

$Q\
\mu^{(\ell,r)}_{t}\\{v:({\mathbb{E}\,}(R^{\gamma})\,t)^{-1/\gamma}v\geq\rho\\}=Q\
\frac{\rho^{-\gamma}}{\gamma}\ (1+o(1))$ (15)

combined with the uniform bound

$\sup_{t>0}\sup_{\rho>0}\rho^{\gamma}\
\mu^{(\ell,r)}_{t}\\{v:({\mathbb{E}\,}(R^{\gamma})\,t)^{-1/\gamma}v\geq\rho\\}<\infty.$
(16)

Indeed, by substituting $\varrho=\rho({\mathbb{E}\,}(R^{\gamma})t)^{1/\gamma}$
in (14) we obtain (15) and by making the same substitution in (13) we obtain
(16). $\Box$

###  3.3 A decomposition

Take some $v_{0}>0$ and split the integral representation (8) into three
parts:

$\displaystyle Y_{Q,\gamma}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{v_{0}}v{\widetilde{N}}_{Q,\gamma}(dv)+\int_{v_{0}}^{\infty}vN_{Q,\gamma}(dv)-Q\int_{v_{0}}^{\infty}v\mu^{(\ell,r)}_{t}(dv)$
(17) $\displaystyle:=$ $\displaystyle Y^{\circ}(t)+Y^{{\dagger}}(t)-E_{t},$

where $N$ is the corresponding non-centered Poisson random measure and $E_{t}$
is the centering deterministic function.

The variance of $Y^{\circ}(t)$ admits an upper bound

$\displaystyle\textrm{Var}\,Y^{\circ}(t)$ $\displaystyle=$ $\displaystyle
Q\int_{0}^{v_{0}}v^{2}\mu^{(\ell,r)}_{t}(dv)=2Q\int_{0}^{v_{0}}v\
\mu^{(\ell,r)}_{t}[v,v_{0}]\,dv$ $\displaystyle\leq$ $\displaystyle
2Q\int_{0}^{v_{0}}v\ \mu^{(\ell,r)}_{t}[v,\infty)\,dv.$

Using (13) we get

$\textrm{Var}\,Y^{\circ}(t)\leq\frac{2\,Q\,t}{\gamma(\gamma-1)}{\mathbb{E}\,}(R^{\gamma})\int_{0}^{v_{0}}v^{1-\gamma}\,dv=D_{2}\,t\,v_{0}^{2-\gamma},$
(18)

where
$D_{2}:=\tfrac{2Q}{\gamma(\gamma-1)(2-\gamma)}\,{\mathbb{E}\,}(R^{\gamma})$.

Similarly, the centering term admits a bound

$0\leq E_{t}\leq
Q\int_{v_{0}}^{\infty}\mu^{(\ell,r)}_{t}[v,\infty)\,dv+Q\,v_{0}\,\mu^{(\ell,r)}_{t}[v_{0},\infty)\leq
D_{1}\,t\,v_{0}^{1-\gamma},$ (19)

where $D_{1}:=\tfrac{Q}{(\gamma-1)^{2}}\,{\mathbb{E}\,}(R^{\gamma})$.

###  3.4 A lower bound for large deviations

We will give a lower bound for large deviation probabilities
${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)$ with $\varrho=\varrho(t)\gg
t^{1/\gamma}$. Let $h,\delta$ be small positive numbers. Define
$v_{0}:=h\varrho$ and consider the corresponding decomposition (17).

First of all, notice that $E_{t}$ is negligible at the range $\varrho$ because
by (19) we have

$E_{t}\leq
D_{1}t(h\varrho)^{1-\gamma}=D_{1}h^{1-\gamma}\left(t^{-1/\gamma}\varrho\right)^{-\gamma}\varrho=o(\varrho).$

Therefore, we may and do assume $t$ to be so large that
$E_{t}\leq\delta\varrho$.

Using (18), by Chebyshev inequality we have

$\displaystyle{\mathbb{P}}(|Y^{\circ}(t)|\geq\delta\varrho)$
$\displaystyle\leq$
$\displaystyle\frac{\textrm{Var}\,Y^{\circ}(t)}{(\delta\varrho)^{2}}\leq\frac{D_{2}t(h\varrho)^{2-\gamma}}{(\delta\varrho)^{2}}$
(20) $\displaystyle=$ $\displaystyle\frac{D_{2}h^{2-\gamma}}{\delta^{2}}\
(t^{-1/\gamma}\varrho)^{-\gamma}\to 0.$

It is also useful to notice that for each $\rho>0$ and all large $t$

$\displaystyle\mu^{(\ell,r)}_{t}[v_{0},\infty)$ $\displaystyle=$
$\displaystyle\mu^{(\ell,r)}_{t}[h\varrho,\infty)=\mu^{(\ell,r)}_{t}[h(t^{-1/\gamma}\varrho)t^{1/\gamma},\infty)$
$\displaystyle\leq$ $\displaystyle\mu^{(\ell,r)}_{t}[h\rho
t^{1/\gamma},\infty)\leq\frac{{\mathbb{E}\,}(R^{\gamma})}{\gamma(\gamma-1)}\,(h\rho)^{-\gamma},$

where we used (13) at the last step. Letting $\rho\to\infty$ we get
$\mu^{(\ell,r)}_{t}[v_{0},\infty)\to 0$.

Now we may proceed with the required lower bound as follows:

$\displaystyle{\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)\geq{\mathbb{P}}(|Y^{\circ}(t)|\leq\delta\varrho,Y^{{\dagger}}(t)\geq(1+2\delta)\varrho)$
$\displaystyle\geq$
$\displaystyle{\mathbb{P}}(|Y^{\circ}(t)|\leq\delta\varrho)\
{\mathbb{P}}(Y^{{\dagger}}(t)\geq(1+2\delta)\varrho;N(v_{0},\infty)=1)$
$\displaystyle=$
$\displaystyle{\mathbb{P}}(|Y^{\circ}(t)|\leq\delta\varrho)\,\exp\\{-Q\mu^{(\ell,r)}_{t}[v_{0},\infty)\\}\,Q\,\mu^{(\ell,r)}_{t}[(1+2\delta)\varrho,\infty).$

The idea behind this bound is to take a single service process providing a
substantial large deviation workload and to suppress other contributions.

As we have just seen, the first two factors tend to one, thus

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)\geq
Q\,\mu^{(\ell,r)}_{t}[(1+2\delta)\varrho,\infty)\ (1+o(1)).$ (21)

###  3.5 An upper bound for large deviations

Starting again with representation (17), using $E_{t}\geq 0$ and (20) we have

$\displaystyle{\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)$ (22)
$\displaystyle\leq$
$\displaystyle{\mathbb{P}}(Y^{\circ}(t)\geq\delta\varrho)+{\mathbb{P}}(N[v_{0},\infty)\geq
2)+{\mathbb{P}}(Y^{{\dagger}}(t)\geq(1-\delta)\varrho;N[v_{0},\infty)=1)$
$\displaystyle=$
$\displaystyle{\mathbb{P}}(Y^{\circ}(t)\geq\delta\varrho)+{\mathbb{P}}(N[v_{0},\infty)\geq
2)+{\mathbb{P}}(N[(1-\delta)\varrho,\infty)=1)$ $\displaystyle\leq$
$\displaystyle\frac{D_{2}th^{2-\gamma}}{\delta^{2}\varrho^{\gamma}}+\frac{1}{2}\left(Q\mu^{(\ell,r)}_{t}[v_{0},\infty)\right)^{2}+Q\mu^{(\ell,r)}_{t}[(1-\delta)\varrho,\infty).$

Here the last term is the main one. Recall that almost the same expression
also shows up in the lower bound.

###  3.6 Proof of Theorem 3

Recall that, according to (14), in the zone under consideration
$t^{1/\gamma}\ll\varrho\ll t$, it is true that

$\mu^{(\ell,r)}_{t}[\varrho,\infty)=\frac{{\mathbb{E}\,}(R^{\gamma})}{\gamma}\
t\ \varrho^{-\gamma}(1+o(1))$ (23)

and we have the similar representations with $\varrho$ replaced by either
$(1+2\delta)\varrho$, $(1-\delta)\varrho$, or $v_{0}=h\varrho$.

In view of (23), the lower estimate (21) yields

$\liminf_{t\to\infty}\frac{{\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)}{t\varrho^{-\gamma}}\geq\frac{Q\
{\mathbb{E}\,}(R^{\gamma})}{\gamma}(1+2\delta)^{-\gamma},$

while the upper estimate (22) yields

$\limsup_{t\to\infty}\frac{{\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)}{t\varrho^{-\gamma}}\leq\frac{D_{2}h^{2-\gamma}}{\delta^{2}}+\frac{Q\
{\mathbb{E}\,}(R^{\gamma})}{\gamma}(1-\delta)^{-\gamma},$

because the second term in (22) has a lower order of magnitude.

Letting first $h\to 0$, then $\delta\to 0$, we obtain

$\lim_{t\to\infty}\frac{{\mathbb{P}}(Y_{Q,\gamma}(t)\geq\varrho)}{t\varrho^{-\gamma}}=\frac{Q\
{\mathbb{E}\,}(R^{\gamma})}{\gamma},$

as required. $\Box$

###  3.7 Proof of Theorem 4

The proof goes along the same lines as in the moderate deviation case, except
for the evaluation of $\mu^{(\ell,r)}_{t}[\varrho,\infty)$. Instead of (23),
we have the following non-asymptotic exact formula. According to (9), for
$\rho=\kappa t$ we have

$\displaystyle\mu^{(\ell,r)}_{t}[\varrho,\infty)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mu^{(\ell)}_{t}\left[\frac{\varrho}{r},\infty\right)F_{R}(dr)$
$\displaystyle=$ $\displaystyle\int_{\kappa}^{\infty}\left(t\
\frac{(\varrho/r)^{-\gamma}}{\gamma}+\frac{2-\gamma}{(\gamma-1)\gamma}\,(\varrho/r)^{1-\gamma}\right)F_{R}(dr)$
$\displaystyle=$
$\displaystyle\int_{\kappa}^{\infty}\left(\frac{\kappa^{-\gamma}}{\gamma}\
r^{\gamma}+\frac{(2-\gamma)\kappa^{1-\gamma}}{(\gamma-1)\gamma}\,r^{\gamma-1}\right)F_{R}(dr)\,t^{-(\gamma-1)}$
$\displaystyle=$
$\displaystyle\left(\frac{\kappa^{-\gamma}}{\gamma}\,{\mathbb{E}\,}(R^{\gamma}{\mathbf{1}}_{\\{R\geq\kappa\\}})+\frac{(2-\gamma)\kappa^{1-\gamma}}{(\gamma-1)\gamma}\,{\mathbb{E}\,}(R^{\gamma-1}{\mathbf{1}}_{\\{R\geq\kappa\\}})\right)\,t^{-(\gamma-1)}$
$\displaystyle=$ $\displaystyle D^{(1)}_{I}(\kappa)\ t^{-(\gamma-1)}.$

The latter constant is positive due to assumption
${\mathbb{P}}(R\geq\kappa)>0$.

For the lower bound, the estimate (21) yields

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\kappa t)\geq
Q\,D^{(1)}_{I}((1+2\delta)\kappa)\,t^{-(\gamma-1)}\ (1+o(1)).$

Letting $\delta\searrow 0$ and using (5), we have

$\displaystyle\lim_{\delta\searrow 0}D^{(1)}_{I}((1+2\delta)\kappa)$
$\displaystyle=$
$\displaystyle\frac{\kappa^{-\gamma}}{\gamma}\,{\mathbb{E}\,}(R^{\gamma}{\mathbf{1}}_{\\{R>\kappa\\}})+\frac{(2-\gamma)\kappa^{1-\gamma}}{(\gamma-1)\gamma}\,{\mathbb{E}\,}(R^{\gamma-1}{\mathbf{1}}_{\\{R>\kappa\\}})$
$\displaystyle=$ $\displaystyle D^{(1)}_{I}(\kappa).$

Therefore,

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\kappa t)\geq
Q\,D^{(1)}_{I}(\kappa)\,t^{-(\gamma-1)}\ (1+o(1)),$

as required.

For the upper bound, the estimate (22) with $\varrho=\kappa t$ yields

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\kappa
t)\leq\left(\frac{D_{2}h^{2-\gamma}}{\delta^{2}\kappa^{\gamma}}+Q\,D^{(1)}_{I}((1-\delta)\kappa)\right)t^{-(\gamma-1)}\
(1+o(1)).$

Letting first $h\searrow 0$, we get rid of the first term and obtain

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\kappa t)\leq
Q\,D^{(1)}_{I}((1-\delta)\kappa)\,t^{-(\gamma-1)}\ (1+o(1)).$

Letting $\delta\searrow 0$, we have

$\lim_{\delta\searrow 0}D^{(1)}_{I}((1-\delta)\kappa)=D^{(1)}_{I}(\kappa).$

Therefore,

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\kappa t)\leq
Q\,D^{(1)}_{I}(\kappa)\,t^{-(\gamma-1)}\ (1+o(1)),$

as required. $\Box$

###  3.8 Proof of Theorem 7

In the setting of this theorem, the large deviation probabilities decay faster
with $t$ than Chebyshev inequality (20) suggests. Therefore, we need a finer
estimate for $Y^{\circ}(t)$ given in the following lemma.

###### Lemma 11

For every $m,\delta>0$ there exist $h>0$ and $C=C(h,\delta)>0$ such that for
all $t>0$

${\mathbb{P}}(Y^{\circ}(t)\geq\delta t)\leq C\,t^{-(\gamma-1)m},$

where $Y^{\circ}(t)$ is defined by $\eqref{split}$ with the splitting point
$v_{0}:=ht$.

Proof of the Lemma: We start with some calculations valid for arbitrary
$v_{0}$. We have the following formula for exponential moment of the centered
Poisson integral:

${\mathbb{E}\,}\exp(\lambda
Y^{\circ}(t))=\exp\big{\\{}\int_{0}^{v_{0}}(e^{\lambda v}-1-\lambda
v)\mu^{(\ell,r)}_{t}(dv)\big{\\}}.$ (24)

Let us split the integration domain in (24) into two parts: $[0,v_{0}/2]$ and
$(v_{0}/2,v_{0}]$. For the second one we have

$\int_{\frac{v_{0}}{2}}^{v_{0}}(e^{\lambda v}-1-\lambda
v)\mu^{(\ell,r)}_{t}(dv)\leq e^{\lambda
v_{0}}\cdot\mu^{(\ell,r)}_{t}[\frac{v_{0}}{2},v_{0}]\leq D_{3}e^{\lambda
v_{0}}tv_{0}^{-\gamma},$ (25)

where
$D_{3}:=\frac{2^{\gamma}\,{\mathbb{E}\,}(R^{\gamma})}{\gamma(\gamma-1)}$. At
the last step we used (13).

For the first zone, by using inequality $e^{x}-1-x\leq x^{2}e^{x}$ and (13) we
have

$\displaystyle\int_{0}^{\frac{v_{0}}{2}}(e^{\lambda v}-1-\lambda
v)\mu^{(\ell,r)}_{t}(dv)\leq\int_{0}^{\frac{v_{0}}{2}}\lambda^{2}v^{2}e^{\lambda
v}\mu^{(\ell,r)}_{t}(dv)$ $\displaystyle\leq$ $\displaystyle 2\,e^{\lambda
v_{0}/2}\lambda^{2}\int_{0}^{\frac{v_{0}}{2}}v\mu^{(\ell,r)}_{t}[v,v_{0}/2]dv\leq\frac{2\,{\mathbb{E}\,}(R^{\gamma})}{\gamma(\gamma-1)}\,e^{\lambda
v_{0}/2}\lambda^{2}t\int_{0}^{\frac{v_{0}}{2}}v^{1-\gamma}dv$ $\displaystyle=$
$\displaystyle D_{4}\,e^{\lambda v_{0}/2}\,\lambda^{2}\,t\,v_{0}^{2-\gamma},$

where
$D_{4}:=\tfrac{2^{\gamma-1}\,{\mathbb{E}\,}(R^{\gamma})}{\gamma(\gamma-1)(2-\gamma)}$.

Next, using inequality $e^{\frac{x}{2}}x^{2}\leq 3e^{x}$, we have

$\int_{0}^{\frac{v_{0}}{2}}(e^{\lambda v}-1-\lambda
v)\mu^{(\ell,r)}_{t}(dv)\leq 3\,D_{4}\,e^{\lambda v_{0}}\,t\,v_{0}^{-\gamma}.$
(26)

By summing up (25) and (26), we obtain

${\mathbb{E}\,}\exp(\lambda
Y^{\circ}(t))\leq\exp\left\\{\left(D_{3}+3\,D_{4}\right)t\,v_{0}^{-\gamma}e^{\lambda
v_{0}}\right\\}:=\exp\left\\{Ae^{\lambda v_{0}}\right\\},$

where

$A:=(D_{3}+3D_{4})tv_{0}^{-\gamma}.$

For every real $y$ by exponential Chebyshev inequality we have

${\mathbb{P}}(Y^{\circ}(t)\geq y)\leq\inf\limits_{\lambda>0}\ \exp(Ae^{\lambda
v_{0}}-\lambda y).$ (27)

If $y>Av_{0}$, the minimum on the right hand side is attained at the point
$\lambda=\frac{1}{v_{0}}\log(\frac{y}{Av_{0}})$. By plugging this value in
(27) we obtain

$\displaystyle{\mathbb{P}}(Y^{\circ}(t)\geq y)$ $\displaystyle\leq$
$\displaystyle\exp\left(\frac{y}{v_{0}}\right)\
\left(\frac{Av_{0}}{y}\right)^{\frac{y}{v_{0}}}$ (28) $\displaystyle=$
$\displaystyle\exp\left(\frac{y}{v_{0}}\right)\
\left((D_{3}+3D_{4})\,\frac{tv_{0}^{1-\gamma}}{y}\right)^{\frac{y}{v_{0}}}.$

Letting here $y:=\delta t$, $v_{0}:=ht$ yields

${\mathbb{P}}(Y^{\circ}(t)\geq\delta t)\leq
C\,t^{-\frac{\delta}{h}(\gamma-1)},$

where $C$ depends only on $\delta,h$. Choosing $h<\frac{\delta}{m}$ we get the
result. $\square$

Now we can proceed to the proof of the theorem.

Upper bound. Let $\eta:=\tfrac{(1-\zeta)\kappa}{n-\zeta}$. Since
$\zeta\in(0,1)$, we have $\eta>0$. It also follows from the definition that
$\tfrac{\kappa}{n-\zeta}=\tfrac{\kappa-\eta}{n-1}$. Therefore, we may rewrite
(6) as

${\mathbb{P}}\left(R\geq\frac{\kappa-\eta}{n-1}\right)=0.$ (29)

Let $\delta\in(0,\eta)$. By using Lemma 11 with $m=n+1$ we find a small $h>0$
such that

${\mathbb{P}}(Y^{\circ}(t)\geq\delta t)\leq C\,t^{-(\gamma-1)(n+1)}.$ (30)

By using the decomposition (17) with $v_{0}=ht$ and taking into account
$E_{t}\geq 0$ we get the bound

${\mathbb{P}}\left(Y_{Q,\gamma}(t)\geq\kappa
t\right)\leq{\mathbb{P}}\left(Y^{\circ}(t)\geq\delta
t\right)+{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(\kappa-\delta)t\right).$

The first term is negligible by (30). Let denote $N_{0}:=N[v_{0},\infty)$,
which is a Poissonian random variable with intensity
$\mu_{0}:=Q\mu^{(\ell,r)}_{t}[v_{0},\infty)$, and apply the following bound to
the second term:

$\displaystyle{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(\kappa-\delta)t\right)\leq{\mathbb{P}}(N_{0}>n)$
$\displaystyle+\
{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(\kappa-\delta)t\right);N_{0}=n)+{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(\kappa-\delta)t\right),N_{0}\leq
n-1).$

For the first term, an elementary bound for Poisson tail works, namely

${\mathbb{P}}(N_{0}>n)=e^{-\mu_{0}}\sum_{m=0}^{\infty}\frac{\mu_{0}^{n+1+m}}{(n+1+m)!}\leq
e^{-\mu_{0}}\frac{\mu_{0}^{n+1}}{(n+1)!}\sum_{m=0}^{\infty}\frac{\mu_{0}^{m}}{m!}\leq\frac{\mu_{0}^{n+1}}{(n+1)!}$

where we used that $(n+1+m)!\geq(n+1)!m!$. Notice that by (13) with
$\rho:=v_{0}=ht$ we have

$\mu_{0}\leq\frac{Q\,{\mathbb{E}\,}(R^{\gamma})}{\gamma(\gamma-1)}\,t\,(ht)^{-\gamma}=\frac{Q\,{\mathbb{E}\,}(R^{\gamma})h^{-\gamma}}{\gamma(\gamma-1)}\,t^{-(\gamma-1)},$

hence,

${\mathbb{P}}(N_{0}>n)=O\left(t^{-(\gamma-1)(n+1)}\right)$

is negligible.

Further, by using (29) and the definition of the measure $\mu^{(\ell,r)}_{t}$
we see that

$\mu^{(\ell,r)}_{t}[\tfrac{(\kappa-\delta)t}{n-1},\infty)\leq\mu^{(\ell,r)}_{t}[\tfrac{(\kappa-\eta)t}{n-1},\infty)=0,$

which implies

${\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(\kappa-\delta)t\right),N_{0}\leq
n-1)=0,$

because here Poissonian integral $Y^{{\dagger}}(t)$ is a sum of not more than
$n-1$ terms each being strictly smaller than $\tfrac{(\kappa-\delta)t}{n-1}$.

For $A\in\mathcal{B}([v_{0},\infty))$ denote $N_{A}:=N(A)$ with intensity
$\mu_{A}:=Q\mu^{(\ell,r)}_{t}(A)$ and
$\nu_{t}^{(l,r)}(A):={\mathbb{P}}(N_{A}=1\ |\ N_{0}=1)$, which is a measure on
$[v_{0},\infty)$. We have

$\nu_{t}^{(l,r)}(A)=e^{-\mu_{A}}\mu_{A}\cdot
e^{\mu_{A}-\mu_{0}}\cdot\frac{e^{\mu_{0}}}{\mu_{0}}=\frac{\mu_{A}}{\mu_{0}}.$

The remaining Poissonian integral with fixed number of points admits the
following representation

$\displaystyle{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(\kappa-\delta)t\right);N_{0}=n)$
$\displaystyle=$
$\displaystyle{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(\kappa-\delta)t\right)\
|\ N_{0}=n)\ {\mathbb{P}}(N_{0})=n$ $\displaystyle=$ $\displaystyle
e^{-\mu_{0}}\
\frac{\mu_{0}^{n}}{n!}\int_{[v_{0},\infty)^{n}}{\mathbf{1}}_{\\{v_{1}+...+v_{n}\geq(\kappa-\delta)t\\}}\prod_{m=1}^{n}\nu_{t}^{(l,r)}(dv_{m})$
$\displaystyle\leq$ $\displaystyle e^{-\mu_{0}}\
\frac{Q^{n}}{n!}\int_{{\mathbb{R}}_{+}^{n}}{\mathbf{1}}_{\\{v_{1}+...+v_{n}\geq(\kappa-\delta)t\\}}\prod_{m=1}^{n}\mu^{(\ell,r)}_{t}(dv_{m})$
$\displaystyle=$ $\displaystyle e^{-\mu_{0}}\
\frac{Q^{n}}{n!}\int_{[0,t]^{n}}\int_{R_{+}^{n}}{\mathbf{1}}_{\\{\ell_{1}r_{1}+...+\ell_{n}r_{n}\geq(\kappa-\delta)t\\}}\prod_{m=1}^{n}\mu^{(\ell)}_{t}(d\ell_{m})\prod_{m=1}^{n}F_{R}(dr_{m})$
$\displaystyle=$ $\displaystyle e^{-\mu_{0}}\
\frac{Q^{n}}{n!}\int_{[0,1]^{n}}{\mathbb{P}}\left(s_{1}R_{1}+...+s_{n}R_{n}\geq\kappa-\delta\right)\prod_{m=1}^{n}\nu(ds_{m})\
t^{-(\gamma-1)n},$ $\displaystyle:=$ $\displaystyle e^{-\mu_{0}}\
Q^{n}D_{I}^{(n)}(\kappa-\delta)\ t^{-(\gamma-1)n},$

where $R_{1},...,R_{n}$ are i.i.d. variables with distribution $F_{R}$ and,
according to (9), $\nu$ is a measure on $[0,1]$ having the weight
$\tfrac{1}{\gamma(\gamma-1)}$ at $1$ and the density

$\frac{d\nu}{ds}\,(s)=s^{-(\gamma+1)}+\frac{2-\gamma}{\gamma}\,s^{-\gamma},\qquad
0<s<1.$

Notice also that the constant $D^{(n)}_{I}(\kappa-\delta)$ is finite although
measure $\nu$ is infinite at each neighborhood of zero. The reason is that the
probability we integrate vanishes if for some $m$ one has
$s_{m}<s_{*}:=\tfrac{(n-1)(\eta-\delta)}{\kappa-\eta}$ where $\eta>0$
satisfies (29). Indeed, we have in this case

$\displaystyle{\mathbb{P}}\left(s_{1}R_{1}+...+s_{n}R_{n}\geq\kappa-\delta\right)$
$\displaystyle\leq$ $\displaystyle{\mathbb{P}}\left((s_{*}+(n-1))\max_{1\leq
m\leq n}R_{j}\geq\kappa-\delta\right)$ $\displaystyle\leq$ $\displaystyle
n\,{\mathbb{P}}\left(R\geq\frac{\kappa-\delta}{s_{*}+(n-1)}\right)$
$\displaystyle=$ $\displaystyle
n\,{\mathbb{P}}\left(R\geq\frac{\kappa-\eta}{n-1}\right)=0.$

We summarize our findings as

${\mathbb{P}}\left(Y_{Q,\gamma}(t)\geq\kappa t\right)\leq
Q^{n}D_{I}^{(n)}(\kappa-\delta)\ t^{-(\gamma-1)n}(1+o(1)).$

Letting $\delta\searrow 0$, we obtain

${\mathbb{P}}\left(Y_{Q,\gamma}(t)\geq\kappa t\right)\leq
Q^{n}\,D^{(n)}_{I}(\kappa)\ t^{-(\gamma-1)n}(1+o(1)),$

where

$\displaystyle D^{(n)}_{I}(\kappa)$ $\displaystyle:=$
$\displaystyle\lim_{\delta\to 0}D^{(n)}_{I}(\kappa-\delta)$ (31)
$\displaystyle=$
$\displaystyle\frac{1}{n!}\int_{[0,1]^{n}}{\mathbb{P}}\left(s_{1}R_{1}+...+s_{n}R_{n}\geq\kappa\right)\prod_{m=1}^{n}\nu(ds_{m}).$

It is easy to see that for $n=1$ we obtain the same value of
$D^{(1)}_{I}(\kappa)$ as in Theorem 4.

Lower bound.

First, notice that $E_{t}$ in (17) is still negligible because by (19) we have

$E_{t}\leq
D_{1}\,t\,v_{0}^{1-\gamma}=D_{1}\,t\,(ht)^{1-\gamma}=O(t^{2-\gamma})=o(t).$

Hence, for every fixed small $\delta$ we may and do assume that
$E_{t}\leq\delta t$ for large $t$.

Second, using (18), by Chebyshev inequality we have

$\displaystyle{\mathbb{P}}(|Y^{\circ}(t)|\geq\delta\,t)$ $\displaystyle\leq$
$\displaystyle\frac{\textrm{Var}\,Y^{\circ}(t)}{(\delta\,t)^{2}}\leq\frac{D_{2}t(ht)^{2-\gamma}}{(\delta\,t)^{2}}$
$\displaystyle=$ $\displaystyle\frac{D_{2}h^{2-\gamma}}{\delta^{2}}\
t^{-(\gamma-1)}\to 0,\quad\textrm{as }t\to\infty.$

Therefore, we may proceed towards the required lower bound as follows:

$\displaystyle{\mathbb{P}}(Y_{Q,\gamma}(t)\geq\kappa\,t)\geq{\mathbb{P}}(|Y^{\circ}(t)|\leq\delta\,t,Y^{{\dagger}}(t)\geq(\kappa+2\delta)t)$
$\displaystyle\geq$ $\displaystyle{\mathbb{P}}(|Y^{\circ}(t)|\leq\delta\,t)\
{\mathbb{P}}(Y^{{\dagger}}(t)\geq(\kappa+2\delta)t;N_{0}=n)$ $\displaystyle=$
$\displaystyle(1+o(1))\
{\mathbb{P}}(Y^{{\dagger}}(t)\geq(\kappa+2\delta)t;N_{0}=n).$

The idea behind this bound is to focus on $n$ service process providing a
substantial large deviation workload and to suppress other contributions.

Furthermore, by using the expression obtained while working on the upper bound

$\displaystyle{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(\kappa+2\delta)t\right);N_{0}=n)$
$\displaystyle=$
$\displaystyle\frac{Q^{n}}{n!}\int_{[0,1]^{n}}{\mathbb{P}}\left(s_{1}R_{1}+...+s_{n}R_{n}\geq\kappa+2\delta\right)\prod_{m=1}^{n}\nu(ds_{m})\
t^{-(\gamma-1)n}(1+o(1))$ $\displaystyle:=$ $\displaystyle
Q^{n}\,D_{I}^{(n)}(\kappa+2\delta)t^{-(\gamma-1)n}(1+o(1)).$

By letting $\delta\searrow 0$, we obtain

$\displaystyle\lim_{\delta\searrow 0}D_{I}^{(n)}(\kappa+2\delta)$
$\displaystyle=$
$\displaystyle\frac{1}{n!}\int_{[0,1]^{n}}{\mathbb{P}}\left(s_{1}R_{1}+...+s_{n}R_{n}>\kappa\right)\prod_{m=1}^{n}\nu(ds_{m})$
$\displaystyle=$
$\displaystyle\frac{1}{n!}\int_{[0,1]^{n}}{\mathbb{P}}\left(s_{1}R_{1}+...+s_{n}R_{n}\geq\kappa\right)\prod_{m=1}^{n}\nu(ds_{m})$
$\displaystyle=$ $\displaystyle D_{I}^{(n)}(\kappa).$

For the non-obvious passage we have used the following lemma.

###### Lemma 12

Assume that (7) holds. Then

$\nu^{n}\left\\{{\mathbf{s}}:=(s_{1},...,s_{n}):{\mathbb{P}}(s_{1}R_{1}+...+s_{n}R_{n}=\kappa)>0\right\\}=0.$
(32)

The required lower bound

${\mathbb{P}}(Y_{Q,\gamma}(t)\geq\kappa\,t)\geq
Q^{n}D^{(n)}_{I}(\kappa)\,t^{-(\gamma-1)n}(1+o(1))$

follows now from the previous estimates. It merely remains to prove the lemma.
$\Box$

Proof of Lemma 12 : Let $r_{1},...,r_{n}$ be the atoms of the distribution
$F_{R}$, i.e. ${\mathbb{P}}(r_{m})>0,1\leq m\leq n$. Define

$F=F(r_{1},...,r_{n}):=\\{{\mathbf{s}}\in[0,1]^{n}:s_{1}r_{1}+...+s_{n}r_{n}=\kappa\\}.$

For every subset of integers $M\subset[1..n]$ let

$B_{M}:=\\{{\mathbf{s}}\in[0,1]^{n}:s_{m}\in[0,1),m\in M;s_{m}=1,m\not\in
M\\}.$

Notice that $[0,1]^{n}=\cup_{M}B_{M}$. Let

$F_{M}:=F\bigcap B_{m}=\\{{\mathbf{s}}\in B_{M}:\sum_{m\in
M}s_{m}r_{m}=\kappa-\sum_{m\not\in M}r_{m}\\}.$

If $M$ is not empty, then $\nu^{n}(F_{M})=0$ because $\nu$ is absolutely
continuous on $[0,1)$.

If $M$ is empty, then $B_{M}=\\{(1,...,1)\\}$ is a singleton and
$F_{M}=\emptyset$ because otherwise $\sum_{m=1}^{n}r_{m}=\kappa$ which would
contradict to (7).

We conclude that

$\nu^{n}\left(F(r_{1},...,r_{n})\right)=\sum_{M}\nu^{n}(F_{M})=0.$

Since

$\left\\{{\mathbf{s}}:{\mathbb{P}}(s_{1}R_{1}+...+s_{n}R_{n}=\kappa)>0\right\\}\subset\bigcup_{r_{1},...,r_{n}}F(r_{1},...,r_{n})$

and the union is countable, we obtain (32).

$\square$

###  3.9 Proof of Theorem 10

Upper bound.

We take a small $\delta>0$, use decomposition (17) with $v_{0}:=h\varrho$ (a
small $h=h(\delta)$ will be specified later on), and start with a usual bound

$\displaystyle{\mathbb{P}}\left(Y_{Q,\gamma}(t)\geq\varrho\right)\leq{\mathbb{P}}\left(Y^{\circ}(t)\geq\delta\varrho\right)+{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(1-\delta)\varrho\right)$
$\displaystyle\leq$
$\displaystyle{\mathbb{P}}\left(Y^{\circ}(t)\geq\delta\varrho\right)+{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(1-\delta)\varrho;N_{0}=1\right)+P(N_{0}\geq
2).$

To show that the first term is negligible we use estimate (28) with
$y:=\delta\varrho$, $v_{0}:=h\varrho$ and obtain for some $C=C(\delta,h)$

${\mathbb{P}}\left(Y^{\circ}(t)\geq\delta\varrho\right)\leq
C\left(t\varrho^{-\gamma}\right)^{\frac{\delta}{h}}\leq
C\varrho^{-(\gamma-1)\,\frac{\delta}{h}}\ll{\bar{F}_{R}}(\varrho)\leq{\bar{F}_{R}}(\varrho/t)$

whenever $h$ is chosen so small that $(\gamma-1)\,\frac{\delta}{h}>m$.

Subsequent evaluation of $Y^{\dagger}(t)$ requires analysis of the measure
$\mu^{(\ell,r)}_{t}$. By using (11) and (9) we obtain

$\displaystyle\mu^{(\ell,r)}_{t}[v,\infty)$ $\displaystyle=$
$\displaystyle\int_{v/t}^{\infty}\mu^{(\ell)}_{t}\left[\frac{v}{r},t\right]F_{R}(dr)$
$\displaystyle=$
$\displaystyle\int_{v/t}^{\infty}\left(\frac{t\,(r/v)^{\gamma}}{\gamma}+\frac{2-\gamma}{(\gamma-1)\gamma}\,(r/v)^{\gamma-1}\right)F_{R}(dr)$
$\displaystyle=$
$\displaystyle\frac{t\,v^{-\gamma}}{\gamma}\int_{v/t}^{\infty}r^{\gamma}F_{R}(dr)+\frac{(2-\gamma)v^{1-\gamma}}{(\gamma-1)\gamma}\int_{v/t}^{\infty}r^{\gamma-1}F_{R}(dr).$

Since the tail of $F_{R}$ is regularly varying, we have the following
asymptotics for the integrals

$\displaystyle\int_{z}^{\infty}r^{\gamma}F_{R}(dr)$ $\displaystyle=$
$\displaystyle\frac{mz^{\gamma}}{m-\gamma}\,{\bar{F}_{R}}(z)\,(1+o(1)),$
$\displaystyle\int_{z}^{\infty}r^{\gamma-1}F_{R}(dr)$ $\displaystyle=$
$\displaystyle\frac{mz^{\gamma-1}}{m-\gamma+1}\,{\bar{F}_{R}}(z)\,(1+o(1)),\quad\textrm{as
}z\to\infty.$

Therefore, we obtain

$\displaystyle\mu^{(\ell,r)}_{t}[v,\infty)$ $\displaystyle=$ $\displaystyle
t^{-(\gamma-1)}\left[\frac{m}{\gamma(m-\gamma)}+\frac{(2-\gamma)m}{(\gamma-1)\gamma(m-\gamma+1)}\right]{\bar{F}_{R}}(v/t)\,(1+o(1))$
(33) $\displaystyle=$
$\displaystyle\frac{m(m-1)}{\gamma(\gamma-1)(m-\gamma+1)(m-\gamma)}\,t^{-(\gamma-1)}{\bar{F}_{R}}(v/t)\,(1+o(1)).$
$\displaystyle=$ $\displaystyle
D\,t^{-(\gamma-1)}{\bar{F}_{R}}(v/t)\,(1+o(1)),\qquad\textrm{as }v\gg t.$

Now the evaluation of $Y^{\dagger}$ is straightforward. Indeed, by (33)

$\displaystyle{\mathbb{P}}\left(Y^{{\dagger}}(t)\geq(1-\delta)\varrho;N_{0}=1\right)\leq
Q\,\mu^{(\ell,r)}_{t}[(1-\delta)\varrho,\infty)$ $\displaystyle=$
$\displaystyle Q\,D\,t^{\gamma-1}{\bar{F}_{R}}((1-\delta)\varrho/t)\,(1+o(1))$
$\displaystyle=$ $\displaystyle
Q\,D\,t^{\gamma-1}{\bar{F}_{R}}(\varrho/t)(1-\delta)^{-m}\,(1+o(1))$

and

$\displaystyle P(N_{0}\geq 2)\leq
Q^{2}\,\mu^{(\ell,r)}_{t}[h\varrho,\infty)^{2}$ $\displaystyle=$
$\displaystyle Q^{2}\,(Dt^{-(\gamma-1)}{\bar{F}_{R}}(h\varrho/t))^{2}(1+o(1))$
$\displaystyle=$ $\displaystyle
Q^{2}\,(Dt^{-(\gamma-1)}{\bar{F}_{R}}(\varrho/t)h^{-m})^{2}(1+o(1))\ll
t^{-(\gamma-1)}{\bar{F}_{R}}(\varrho/t).$

By combining these estimates and letting $\delta\to 0$ we obtain the desired
bound

${\mathbb{P}}\left(Y_{Q,\gamma}(t)\geq\varrho\right)\leq
Q\,D\,t^{-(\gamma-1)}{\bar{F}_{R}}(\varrho/t)\,(1+o(1)).$

Lower bound.

Since $\varrho\gg t\gg t^{1/\gamma}$, all bounds from section 3.4 apply. For
every $\delta>0$ inequality (21) along with (33) yield

$\displaystyle{\mathbb{P}}\left(Y_{Q,\gamma}(t)\geq\varrho\right)$
$\displaystyle\geq$ $\displaystyle
Q\,\mu^{(\ell,r)}_{t}[(1+2\delta)\varrho,\infty)(1+o(1))$ $\displaystyle=$
$\displaystyle
Q\,D\,t^{-(\gamma-1)}(1+2\delta)^{-m}{\bar{F}_{R}}(\varrho/t)\,(1+o(1)),$

and letting $\delta\to 0$ we get the desired bound

${\mathbb{P}}\left(Y_{Q,\gamma}(t)\geq\varrho\right)\geq
Q\,D\,t^{-(\gamma-1)}{\bar{F}_{R}}(\varrho/t)\,(1+o(1)).$

$\Box$

## Acknowledgement.

This work was supported by Russian Science Foundation grant 21-11-00047.

## References

  * [1]
  * [2] Cohen, S. and Taqqu, M. (2004). Small and large scale behavior of the Poissonized Telecom Process, Methodol. Comput. Appl. Probab. 6, pp. 363–379.
  * [3] Dembo, A., Zeitouni, O. Large Deviations Techniques and Applications, Springer, 2010.
  * [4] Gaigalas, R. (2006). A Poisson bridge between fractional Brownian motion and stable Lévy motion, Stoch. Proc. Appl. 116, pp. 447–462.
  * [5] Kaj, I. (2002). Stochastic Modeling in Broadband Communications Systems, SIAM Monographs on Mathematical Modeling and Computation Vol.8 (SIAM, Philadelphia).
  * [6] Kaj, I. (2005). Limiting fractal random processes in heavy-tailed systems, Fractals in Engineering, New Trends in Theory and Applications, J. Levy-Vehel, E. Lutton (eds.), pp. 199–218 (Springer-Verlag, London).
  * [7] Kaj, I. (2006). Aspects of Wireless Network Modeling Based on Poisson Point Processes, Fields Institute Workshop on Applied Probability (Carleton University, Ottawa).
  * [8] Kaj, I., Leskelä, L., Norros, I., and Schmidt, V. (2007). Scaling limits for random fields with long-range dependence, Ann. Probab. 35, pp. 528–550.
  * [9] Kaj, I. and Taqqu, M. S. (2008). Convergence to fractional Brownian motion and to the Telecom process: the integral representation approach, In and Out of Equilibrium. II., ser.: Progress in Probability, Vol. 60,(Birkhäuser, Basel), pp. 383–427.
  * [10] Kurtz, T. G. (1996). Limit theorems for workload input models, Stochastic Networks, Theory and Applications (Clarendon Press, Oxford) Kelly, F. P., Zachary, S. and Ziedins, I. (eds.), pp. 119–140.
  * [11] Lifshits, M. (2014) Random Processes by Example. World Scientific, Singapore.
  * [12] Pipiras V. and Taqqu, M. S. (2000). The limit of a renewal-reward process with heavy-tailed rewards is not a linear fractional stable motion, Bernoulli 6, pp. 607–614.
  * [13] Rosenkrantz, W. A. and Horowitz, J. (2002). The infinite sourse model for internet traffic: statistical analysis and limit theorems, Methods and Applications of Analysis 9, pp. 445–462.
  * [14] Taqqu, M. S. (2002). The modeling of Ethernet data and of signals that are heavy-tailed with infinite variance, Scand. J. Statist. 29, pp. 273–295.
  * [15]


11institutetext: Holger Götz 22institutetext: Thorsten Pöschel33institutetext:
Institute for Multiscale Simulations
Friedrich-Alexander-Universität Erlangen-Nürnberg
Cauerstraße 3, 91058 Erlangen
Germany
33email<EMAIL_ADDRESS>
# DEM-Simulation of thin elastic membranes interacting with a granulate

Holger Götz  Thorsten Pöschel

(Received: / Accepted: date)

###### Abstract

For a wide range of applications, we need DEM simulations of granular matter
in contact with elastic flexible boundaries. We present a novel method to
describe the interaction between granular particles and a flexible elastic
membrane. Here, the standard mass-spring model approach is supplemented by
surface patches given by a triangulation of the membrane. In contrast to
standard mass-spring models, our simulation method allows for an efficient
simulation even for large particle size dispersion. The novel method allows
coarsening of the mass-spring system leading to a substantial increase of
computation efficiency. The simulation method is demonstrated and benchmarked
for a triaxial test.

###### Keywords:

DEM boundary condition elastic membrane triaxial test

††journal: Granular Matter

##  1 Introduction

For many applications in granular matter research, the system boundaries are
given by deformable containers which may be modeled as elastic membranes. Of
particular interest are jamming systems where the granulate changes its
mechanical properties drastically when the particle number density in the
system is changed liuJammingNotJust1998 ; ohernJammingZeroTemperature2003 ;
ciamarraRecentResultsJamming2010 , which is frequently achieved by evacuating
the air from a deformable container partly filled by granulate. Prominent
examples are granular robotic grippers brownUniversalRoboticGripper2010 ,
where this mechanism is used to grip and manipulate objects, granular paws
hauserFrictionDampingCompliant2016  and similar
fitzgeraldReviewJammingActuation2020 .

In these cases, the system dynamics is determined by two-way coupling, that
is, the deformation of the membrane (e.g. caused by external air pressure)
implies forces on the granular particles and the membrane is deformed under
the action of the granular packing. Such membranes can be modeled as mass-
spring systems (MSS), e.g. debonoDiscreteElementModelling2012 . Here the
membrane is described as a regular or irregular graph whose vertices are
particles and whose edges are linear or non-linear elastic springs. The
elastic behavior is, thus, described by the springs and the contacts between
the membrane and the enclosed granular particles are described by the contacts
between the membrane’s particles and the particles of the granulate.

The choice of the membrane’s structure is critical: If the meshes are too
coarse, particles can penetrate the membrane. Therefore, the mesh width has to
be chosen according to the smallest particles in the system
quDiscreteElementModelling2019 ; debonoDiscreteElementModelling2012 . This is
problematic for several reasons: First, in the coarse of the simulation when
the membrane and the granulate particles interact, the mesh width may change
which is difficult to predict. Second, for the simulation of a highly disperse
granulate, the number of membrane particles and springs can be very large,
resulting in an inefficient simulation. The problem can be solved by
artificially increasing the sizes of the non-interacting membrane particles
(overlapping particles) wuStudyShearBehavior2021  which, however, introduces
an undesired thickness of the membrane. All MSS models are problematic for the
modeling of tangential (frictional) forces along the membrane since the
contacts of the membrane particles - granular particles depends on the
concrete arrangement of the particle positions.

In the current paper, we describe a novel type of MSS which allows for the
simulation of impenetrable flexible elastic boundaries requiring a moderate
number of membrane particles. This model does not suffer from the drawbacks
discussed above. Since by construction the membrane is impenetrable, the
meshes can be chosen larger which makes our method computationally efficient.
The proposed model was used recently to simulate a granular gripper
gotzSoftParticlesReinforce2021  and a bending beam of granular meta-material
bendingBeamMeta:2022 .

##  2 Model description

###  2.1 Particle: particle interaction

The discrete element method (DEM) solves Newton’s equation for the position
$\boldsymbol{r}_{i}$ and the angular orientation $\boldsymbol{\varphi}_{i}$ of
each particle, $i$, of mass $m_{i}$ and tensorial moment of inertia,
$\hat{J}_{i}$:

$\displaystyle m_{i}\frac{\text{d}^{2}\boldsymbol{r}_{i}}{\text{d}t^{2}}$
$\displaystyle=\boldsymbol{F}_{i}=\sum\limits_{j\neq
i}\boldsymbol{F}_{ij}+\boldsymbol{F}_{i}^{\text{ext}}\,,$ (1)
$\displaystyle\hat{J}_{i}\frac{\text{d}^{2}\boldsymbol{\varphi}_{i}}{\text{d}t^{2}}$
$\displaystyle=\boldsymbol{M}_{i}=\sum\limits_{j\neq i}\boldsymbol{M}_{ij}\,.$
(2)

Here, $\boldsymbol{F}_{i}^{\text{ext}}$ is an external force, e.g. gravity,
and $\boldsymbol{F}_{ij}$ and $\boldsymbol{M}_{ij}$ are the force and torque
acting on particle $i$ due to contacts with particles $j$. There are several
models for $\boldsymbol{F}_{ij}$ and $\boldsymbol{M}_{ij}$ as functions of the
relative position, orientation, velocity and angular velocity of the involved
particles, $i$ and $j$, for an extended discussion see, e.g.,
shaferForceSchemesSimulations1996 ; poschelComputationalGranularDynamics2005 ;
kruggel-emdenReviewExtensionNormal2007 ; kruggel-emdenStudyTangentialForce2008
; matuttisUnderstandingDiscreteElement2014 .

###  2.2 Membrane: topology

In the current paper, we focus on the description of an ambient membrane and
its interaction with the granular particles. In our model, the elastically
deformable membrane is modeled by mass-carrying particles that are connected
by viscoelastic springs (mass-spring system, MSS). The topology of the
membrane is given by a mathematical graph whose vertices and edges are
represented by particles and springs, respectively.

The positions of the membrane particles (here called vertex particles), thus,
describe the shape of the membrane (Fig. 1).

Figure 1: Sketch of the MSS

They are subject of Newton’s equations where the forces acting on the vertex
membrane particles originate from three contributions: (a) viscoelastic
stretching of the membrane, (b) moments due to bending of the membrane, and
(c) interaction of granular particles with the membrane. We shall discuss
these contributions in Secs. 2.3-2.5. The total force acting on the vertex
particles is the sum of these three contributions.

###  2.3 Membrane: stretching

Given two adjacent vertex particles $i,j$ at positions
$\boldsymbol{\rho}_{i}$, $\boldsymbol{\rho}_{j}$ and velocities
$\dot{\boldsymbol{\rho}}_{i}$, $\dot{\boldsymbol{\rho}}_{j}$, we define the
relative quantities

$\displaystyle\boldsymbol{\rho}_{ij}$
$\displaystyle\equiv\boldsymbol{\rho}_{i}-\boldsymbol{\rho}_{j},$ (3)
$\displaystyle\dot{\boldsymbol{\rho}}_{ij}$
$\displaystyle\equiv\dot{\boldsymbol{\rho}}_{i}-\dot{\boldsymbol{\rho}}_{j}$
(4)

and the unit vector

$\displaystyle\hat{\boldsymbol{\rho}}_{ij}$
$\displaystyle\equiv\frac{\boldsymbol{\rho}_{ij}}{\rho_{ij}}\,.$ (5)

The interaction between particles $i$ and $j$ is due to linear elastic spring.
Particle $i$ feels the force

$\displaystyle\boldsymbol{F}_{ij}^{\,\text{spring}}=\hat{\boldsymbol{\rho}}_{ij}\left[k\left(\rho_{ij}-\rho_{ij}^{0}\right)-2\gamma\sqrt{k\,m_{ij}^{\text{eff}}}\,\hat{\boldsymbol{\rho}}_{ij}\cdot\dot{\boldsymbol{\rho}}_{ij}\right],$
(6)

which is the force of a damped harmonic oscillator, with equilibrium length
$\rho_{ij}^{0}$, effective mass
$m_{ij}^{\text{eff}}=m_{i}m_{j}/(m_{i}+m_{j})$, damping coefficient $\gamma$,
and spring constant $k$.

To relate the spring constant, $k$, to material characteristics, we notice
that each realistic membrane has a finite width, $d$, and the elasticity of
the membrane material is characterized by its elastic modulus $E$. For our
idealized two-dimensional membrane of vanishing thickness, one obtains the
elastic constant kotMassSpringModels2017 ; ostoja-
starzewskiLatticeModelsMicromechanics2002 .

$k=\frac{\sqrt{3}}{2}E\,d\,.$ (7)

###  2.4 Membrane: flexibility

To explain the description of membrane flexibility, we consider four adjacent
vertex particles at positions $\boldsymbol{\rho}_{1}$,
$\boldsymbol{\rho}_{2}$, $\boldsymbol{\rho}_{3}$, $\boldsymbol{\rho}_{4}$
bridsonSimulationClothingFolds2003 , see Fig. 2.

Figure 2: The deformation of the membrane is described by the angle
$\theta_{12}$ between the normal vectors of adjacent triangles $\triangle 134$
and $\triangle 243$. $\theta_{12}$ is extremely exagerated in this sketch.

The vertex particles span two triangles with normal vectors

$\displaystyle\boldsymbol{n}_{1}$
$\displaystyle=(\boldsymbol{\rho}_{1}-\boldsymbol{\rho}_{3})\times(\boldsymbol{\rho}_{1}-\boldsymbol{\rho}_{4})$
(8) $\displaystyle\boldsymbol{n}_{2}$
$\displaystyle=(\boldsymbol{\rho}_{2}-\boldsymbol{\rho}_{4})\times(\boldsymbol{\rho}_{2}-\boldsymbol{\rho}_{3})\,.$
(9)

The corresponding angle $\theta_{12}$

$\displaystyle\cos\theta_{12}\equiv\boldsymbol{n}_{1}\cdot\boldsymbol{n}_{2}$
(10)

characterizes the flection of the triangles, with respect to their common edge
$\boldsymbol{\rho}_{43}=\boldsymbol{\rho}_{4}-\boldsymbol{\rho}_{3}$. The
restoring torque counteracting the flection can be expressed by elastic and
dissipative forces, $\boldsymbol{F}_{i}^{\text{el}}$, acting on the involved
vertex particles, $i\in\\{1,2,3,4\\}$ bridsonSimulationClothingFolds2003 ,

$\displaystyle\boldsymbol{F}_{i}^{\text{el}}$
$\displaystyle=k^{\text{el}}\frac{\left|\boldsymbol{\rho}_{43}\right|^{2}}{|\boldsymbol{n}_{1}|+|\boldsymbol{n}_{2}|}\left(\sin{\frac{\theta_{12}}{2}}-\sin{\frac{\theta_{12}^{0}}{2}}\right)\boldsymbol{u}_{i},$
(11) $\displaystyle\boldsymbol{F}_{i}^{\text{diss}}$
$\displaystyle=-k^{\text{diss}}\,|\boldsymbol{\rho}_{43}|\,\dot{\theta}_{12}\boldsymbol{u}_{i}\,,$
(12)

where $k^{\text{el}}$ and $k^{\text{diss}}$ are material parameters.

The directions of the forces are linear combinations of the triangles’ normal
vectors given by

$\displaystyle\boldsymbol{u}_{1}$
$\displaystyle=|\boldsymbol{\rho}_{43}|\frac{\boldsymbol{n}_{1}}{|\boldsymbol{n}_{1}|^{2}},$
(13) $\displaystyle\boldsymbol{u}_{2}$
$\displaystyle=|\boldsymbol{\rho}_{43}|\frac{\boldsymbol{n}_{2}}{|\boldsymbol{n}_{2}|^{2}},$
(14) $\displaystyle\boldsymbol{u}_{3}$
$\displaystyle=\frac{(\boldsymbol{x}_{1}-\boldsymbol{x}_{4})\cdot\boldsymbol{\rho}_{43}}{|\boldsymbol{\rho}_{43}|}\frac{\boldsymbol{n}_{1}}{|\boldsymbol{n}_{1}|^{2}}+\frac{(\boldsymbol{x}_{2}-\boldsymbol{x}_{4})\cdot\boldsymbol{\rho}_{43}}{|\boldsymbol{\rho}_{43}|}\frac{\boldsymbol{n}_{2}}{|\boldsymbol{n}_{2}|^{2}},$
(15) $\displaystyle\boldsymbol{u}_{4}$
$\displaystyle=\frac{(\boldsymbol{x}_{1}-\boldsymbol{x}_{3})\cdot\boldsymbol{\rho}_{43}}{|\boldsymbol{\rho}_{43}|}\frac{\boldsymbol{n}_{1}}{|\boldsymbol{n}_{1}|^{2}}+\frac{(\boldsymbol{x}_{2}-\boldsymbol{x}_{3})\cdot\boldsymbol{\rho}_{43}}{|\boldsymbol{\rho}_{43}|}\frac{\boldsymbol{n}_{2}}{|\boldsymbol{n}_{2}|^{2}}\,.$
(16)

Each vertex particle is involved in 6 different pairs of triangles, see Fig.
1. The total force acting on a vertex particle is, thus, the sum of the 6
corresponding forces given by Eqs. (11, 12).

###  2.5 Membrane: granulate-membrane interaction

For the description of the interaction between the membrane and the confined
granular particles, we assume triangular patches spanned between the time
dependent momentary positions of adjacent vertex particles, Fig. 3.

Figure 3: Sketch of a contact between a granular particle and a triangular
patch. The resulting force $\boldsymbol{F}_{c}$ is mapped to the involved
vertex particles $\boldsymbol{\rho}_{1}$, $\boldsymbol{\rho}_{2}$,
$\boldsymbol{\rho}_{3}$ according to the barycentric weights $a_{1}$, $a_{2}$,
$a_{3}$ of the contact point $\boldsymbol{x}_{c}$ with respect to the
locations of the vertex particles.

The interaction of the granular particles with the membrane is then described
by contacts between the granular particles and the patches. This assures that
the patches are always impenetrable disregarding of the sizes of the particles
and the deformation of the membrane.

Contacts between a patch and a granular particle are classified as vertex,
edge, or face contact:

$\begin{split}\text{vertex contact if
~{}}&\boldsymbol{x}_{c}=\boldsymbol{\rho}_{i}\,;~{}i\in\\{1,2,3\\}\\\
\text{edge contact
if~{}}&\left(\boldsymbol{x}_{c}-\boldsymbol{\rho}_{j}\right)\times\left(\boldsymbol{\rho}_{i}-\boldsymbol{\rho}_{j}\right)=\boldsymbol{0}\,;\\\
&~{}i,\,j\in\\{1,2,3\\}\,;~{}i\neq j\\\ \text{face
contact~{}}&\text{else}\end{split}$ (17)

In case a granular particle contacts two neighboring patches, $A$ and $B$, we
chose the adequate contacts according to Tab. 1 in dependence whether these
contacts are face, edge, or vertex contacts huNewAlgorithmContact2013

| face $A$ | edge $A$ | vertex $A$
---|---|---|---
face $B$ |  $A$ and $B$ | $B$ | $B$
edge $B$ | $A$ |  $A$ or $B$ (rand) | $B$
vertex $B$ | $A$ | $A$ |  $A$ or $B$ (rand)
Table 1: Selection of contacts for the case that a granular particle ist in
contact with two adjacent patches, $A$ and $B$

These selection rules do only apply if the patches $A$ and $B$ have a common
edge. Otherwise, all contacts are handles regularly.

In case a granular particle contacts three neighboring patches at their common
vertex, one of these contacts is selected randomly and the others are
disregarded.

Contacts of the membrane with itself may be calculated from contacts between
vertex particles and triangular patches.

Once the contact point is defined, we compute the force according to the
specified contact law. The relative velocity of the granular particle and the
membrane at the contact point which enters the force is interpolated from the
velocities of the vertex particles using barycentric weights, as sketched in
Fig. 3. Similarly, the obtained force is distributed to the involved vertex
particles with barycentric weights $a_{1}$, $a_{2}$, $a_{3}$. The positions
and velocities of the vertex particles and, thus, the dynamics of the membrane
are obtained by numerical integration in the same way as the granular
particles.

The above selection rule above in combination with the barycentric partition
of the force leads to smooth and physically plausible forces acting on the
vertex particles. huNewAlgorithmContact2013

##  3 Applications

We implemented the described flexible wall into the DEM program MercuryDPM
weinhartFastFlexibleParticle2020 . Here we present two examples of its
application.

###  3.1 Triaxial test

The triaxial test is commonly used to investigate the mechanical properties of
a granular sample. To that aim, the sample is placed between two parallel
platens and wrapped by a cylindrical membrane, see Fig. 4.

Figure 4: A triaxial test cell in its initial and final states

A confining stress $\sigma_{\text{c}}$ in radial direction is applied through
the membrane. By controlling the initial distance $h=h_{0}$ between the
platens, an initial stress $\sigma_{\text{p}}=\sigma_{\text{c}}$ is applied in
axial direction. Then, the platens are displaced at constant relative velocity
$\boldsymbol{v_{p}}$ to apply a strain $\epsilon=\log(h_{0}/h)$. We record the
corresponding deviatoric stress
$\sigma_{\text{d}}=\sigma_{\text{p}}-\sigma_{\text{c}}$.

In the simulation, we represent the platens by rigid walls and the membrane by
an MSS. We place small particles at random positions inside the membrane such
that they do not contact one another. To generate the initial conditions, we
apply the Lubachevsky–Stillinger algorithm
lubachevskyGeometricPropertiesRandom1990 . We then gradually apply an initial
stress, $\sigma_{p}$, to the platens and a confining stress,
$\sigma_{\text{c}}$, to the membrane. To this end, each triangular patch of
area $A^{w}_{i}$ and normal unit vector $\hat{\boldsymbol{n}}^{\,w}_{i}$ is
loaded with a force

$\boldsymbol{F}^{\,w}_{i}=A^{\,w}_{i}\sigma_{\text{c}}\,\hat{\boldsymbol{n}}^{\,w}_{i}\,,$
(18)

where $\hat{\boldsymbol{n}}^{\,w}_{i}$ is defined such that force acts from
outside the membrane to the granulate located inside. After this
initialization, we displace the platens at relative velocity
$\boldsymbol{v}_{\text{p}}$ and record the deviatoric stress,
$\sigma_{\text{d}}$. Figure 4 shows the initial and final states of such a
simulation.

To demonstrate the performance of our model, we perform simulations for two
different cases:

  1. 1.

We describe the membrane by a MSS where the confining stress of the membrane
is provided by contacts between vertex particles and granular particles. This
approach was used before, e.g., in quDiscreteElementModelling2019 . We term
this setup _$\text{MSS}_{\text{particle}}$_.

  2. 2.

We describe the membrane as described in Sec. 2. Here the confining stress of
the membrane is provided by contacts between the granular particles and the
triangular patches. We term this setup _$\text{MSS}_{\text{patch}}$_.

For $\text{MSS}_{\text{particle}}$, we place the vertex particles of the
membrane spaced by $s=1.33\,\text{mm}$, such that the radius of a vertex
particle is about $1/3$ the radius of a granular particle. This value is a
trade-off between keeping the computational cost low and having a smooth
membrane that prevents penetration of granular particles
debonoDiscreteElementModelling2012 ; quDiscreteElementModelling2019 . For
$\text{MSS}_{\text{patch}}$, we use $s=4.0\,\text{mm}$ because the membrane
has a smooth surface and is impenetrable by design.

We perform the triaxial test at velocity
$v_{\text{p}}=0.05\,\text{m}/\text{s}$ and confining pressure
$\sigma_{\text{c}}=100\,\text{kPa}$. The material parameters are given in Tab.
2. Furthermore, we choose the parameters $\gamma\approx 0.15$,
$k^{\text{el}}=10^{-3}\,\text{N/m}$ and
$k^{\text{diss}}=10^{-4}\,\text{Ns/m}$. Using the initial areas of the
triangles, we use Eqs. (11-16) to compute the confining forces.

Table 2: Material parameters used in the simulation | particles | membrane
---|---|---
radius / thickness [mm] |  $2.3$ to $2.7$ | 0.3
elastic modulus [Pa] | $4.6\cdot 10^{10}$ | $1.25\cdot 10^{6}$
Poisson’s ratio | 0.245 | 1/3
friction coefficient | 0.25 | 1.2

We consider two different systems: _case 1_ – a cylinder of initial height
$100\,\text{mm}$ and radius $25\,\text{mm}$, and _case 2_ \- a cylinder of
initial height $140\,\text{mm}$ and radius $35\,\text{mm}$. The numbers of
particles used for the membrane and the granulate are given in Tab. 3.

Table 3: Number of granular and membrane particles in the simulations | membrane particles | granular particles
---|---|---
| case 1 | case 2 | case 1 | case 2
$\text{MSS}_{\text{particle}}$ | 10440 | 20374 | $1690$ | $4790$
$\text{MSS}_{\text{patch}}$ | 1131 | 2214 | $1720$ | $4670$

Figure 5

(a) Case 1. initial height $100\,\text{mm}$; radius $25\,\text{mm}$ (b) Case
2. initial height $140\,\text{mm}$; radius $35\,\text{mm}$

Figure 5: Deviatoric stress against axial strain for two different geometries

shows the deviatoric stress, $\sigma_{\text{d}}$, as a function of the axial
strain, $\epsilon$. The different membrane representations,
$\text{MSS}_{\text{particle}}$ and $\text{MSS}_{\text{patch}}$, do not lead to
significant differences of the stress-strain behavior.

Table 4 compares the computer time used to simulate a real time of
$20\,\text{ms}$.

Table 4: Computer time used for different membrane representations. Real time is $20\,\text{ms}$ | case 1 | case 2
---|---|---
$\text{MSS}_{\text{particle}}$ | $\approx 32\,\text{min}$ | $\approx 96\,\text{min}$
$\text{MSS}_{\text{patch}}$ | $\approx 5\,\text{min}$ | $\approx 17\,\text{min}$

The reduced number of particles in $\text{MSS}_{\text{patch}}$ accelerates the
simulations by about the factor 5, compared to simulations using
$\text{MSS}_{\text{particle}}$.

###  3.2 Friction test

A correct representation of frictional forces at contacts with the membrane is
important for many applications. For instance, one of the mechanisms allowing
a granular gripper to grasp an object relies on frictional forces
brownUniversalRoboticGripper2010 . By construction, an object gliding on a
membrane modeled by particles cannot feel a constant (or, at least, smnooth)
force.

We demonstrate the smoothness of a membrane modeled by the here described
model, and the resulting frictional forces by means of a simple sliding test:
In $\text{MSS}_{\text{particle}}$, a membrane is modeled by vertex particles
with spacing $1.33\,\text{mm}$. In $\text{MSS}_{\text{patch}}$, the membrane
is modeled by patches of side length $1.33\,\text{mm}$.

For the test, we place a spherical particle of radius $2\,\text{mm}$ on a
membrane, the free motion of this particle is restricted to the vertical
coordinate, perpendicular to the membrane. Its horizontal motion at constant
velocity, $0.05\,\text{m/s}$, is enforced externally. Figure 6 shows the
particle’s vertical position and the frictional force in sliding direction.

Figure 6: The frictional force in horizontal direction felt a particle when
moving on the surface of a membrane (upper image) and corresponding vertical
position (lower image) vs the horizontal position for both considered membrane
models. Rolling degree of freedom was suppressed

For $\text{MSS}_{\text{particle}}$ we see oscillations in both horizontal
force and vertical position. For $\text{MSS}_{\text{patch}}$, we instead
observe constant vertical position and constant friction force, as expected
when a particle slides on an even plane.

##  4 Conclusions

Previous work using MSS for representing membranes within DEM simulations
describe contacts between the granulate and confining membranes through
contacts between the membrane’s particles and the particles of the granulate.
In this paper, we introduced a membrane model using surface patches. Contacts
between granuate and the membrane are described by contacts between surface
patches constituting the membrane and the particles of the granulate.

The movel model describes a closed surface by design, therefore, the number of
particles in the MSS can be reduced drastically. In a sample simulation
modeling a triaxial test, we obtained an acceleration of the mumerical method
by about a factor five. The comparison of our results with the results using a
traditional membrane description did not reveal significant differences in the
physical behavior. We further demonstrated the new model’s ability to
represent smooth surfaces.

###### Acknowledgements.

We gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) – project number 411517575. The work was supported
by the Interdisciplinary Center for Nanostructured Films (IZNF), the
Competence Unit for Scientific Computing (CSC), and the Interdisciplinary
Center for Functional Particle Systems (FPS) at Friedrich-Alexander
Universität Erlangen-Nürnberg.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * (1) Liu, A. J., Nagel, S. R.: Jamming is not just cool any more.  Nature 396, 21–22 (1998).  DOI 10.1038/23819
  * (2) O’Hern, C. S., Silbert, L. E., Liu, A. J., Nagel, S. R.: Jamming at zero temperature and zero applied stress: The epitome of disorder.  Phys. Rev. E 68, 011306 (2003).  DOI 10.1103/PhysRevE.68.011306
  * (3) Ciamarra, M. P., Nicodemi, M., Coniglio, A.: Recent results on the jamming phase diagram.  Soft Matter 6, 2871–2874 (2010).  DOI 10.1039/B926810C
  * (4) Brown, E., Rodenberg, N., Amend, J., Mozeika, A., Steltz, E., Zakin, M. R., Lipson, H., Jaeger, H. M.: Universal robotic gripper based on the jamming of granular material.  PNAS 107, 18809–18814 (2010).  DOI 10.1073/pnas.1003250107
  * (5) Hauser, S., Eckert, P., Tuleu, A., Ijspeert, A.: Friction and damping of a compliant foot based on granular jamming for legged robots.  In: 2016 6th IEEE International Conference on Biomedical Robotics and Biomechatronics (BioRob), pp. 1160–1165 (2016).  DOI 10.1109/BIOROB.2016.7523788
  * (6) Fitzgerald, S. G., Delaney, G. W., Howard, D.: A Review of jamming actuation in soft robotics.  Actuators 9, 104 (2020).  DOI 10.3390/act9040104
  * (7) de Bono, J., Mcdowell, G., Wanatowski, D.: Discrete element modelling of a flexible membrane for triaxial testing of granular material at high pressures.  Geotech. Lett. 2, 199–203 (2012).  DOI 10.1680/geolett.12.00040
  * (8) Qu, T., Feng, Y. T., Wang, Y., Wang, M.: Discrete element modelling of flexible membrane boundaries for triaxial tests.  Comput. Geotech. 115, 103154 (2019).  DOI 10.1016/j.compgeo.2019.103154
  * (9) Wu, K., Sun, W., Liu, S., Zhang, X.: Study of shear behavior of granular materials by 3D DEM simulation of the triaxial test in the membrane boundary condition.  Adv. Powder Technol. 32, 1145–1156 (2021).  DOI 10.1016/j.apt.2021.02.018
  * (10) Götz, H., Santarossa, A., Sack, A., Pöschel, T., Müller, P.: Soft particles reinforce robotic grippers: Robotic grippers based on granular jamming of soft particles.  Granular Matter 24, 31 (2021).  DOI 10.1007/s10035-021-01193-4
  * (11) Götz, H., Pöschel, T.: Granular meta-material: Viscoelastic response of a bending beam.  Granular Matter sumitted (2022)
  * (12) Schäfer, J., Dippel, S., Wolf, D. E.: Force schemes in simulations of granular materials.  J. Phys. I France 6, 5–20 (1996).  DOI 10.1051/jp1:1996129
  * (13) Pöschel, T., Schwager, T.: Computational Granular Dynamics: Models and Algorithms.  Springer-Verlag, Berlin ; New York (2005)
  * (14) Kruggel-Emden, H., Simsek, E., Rickelt, S., Wirtz, S., Scherer, V.: Review and extension of normal force models for the Discrete Element Method.  Powder Technology 171, 157–173 (2007).  DOI 10.1016/j.powtec.2006.10.004
  * (15) Kruggel-Emden, H., Wirtz, S., Scherer, V.: A study on tangential force laws applicable to the discrete element method (DEM) for materials with viscoelastic or plastic behavior.  Chemical Engineering Science 63(6), 1523–1541 (2008).  DOI 10.1016/j.ces.2007.11.025
  * (16) Matuttis, H. G., Chen, J.: Understanding the Discrete Element Method: Simulation of Non-Spherical Particles for Granular and Multi-body Systems.  Wiley (2014)
  * (17) Kot, M., Nagahashi, H.: Mass spring models with adjustable Poisson’s ratio.  Vis. Comput. 33, 283–291 (2017).  DOI 10.1007/s00371-015-1194-8
  * (18) Ostoja-Starzewski, M.: Lattice models in micromechanics.  Appl. Mech. Rev. 55, 35–60 (2002).  DOI 10.1115/1.1432990
  * (19) Bridson, R., Marino, S., Fedkiw, R.: Simulation of clothing with folds and wrinkles.  In: Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’03, pp. 28–36. Eurographics Association, Goslar, DEU (2003)
  * (20) Hu, L., Hu, G. M., Fang, Z. Q., Zhang, Y.: A new algorithm for contact detection between spherical particle and triangulated mesh boundary in discrete element method simulations.  Int. J. Numer. Methods Eng. 94, 787–804 (2013).  DOI 10.1002/nme.4487
  * (21) Weinhart, T., Orefice, L., Post, M., van Schrojenstein Lantman, M. P., Denissen, I. F. C., Tunuguntla, D. R., Tsang, J. M. F., Cheng, H., Shaheen, M. Y., Shi, H., Rapino, P., Grannonio, E., Losacco, N., Barbosa, J., Jing, L., Alvarez Naranjo, J. E., Roy, S., den Otter, W. K., Thornton, A. R.: Fast, flexible particle simulations — An introduction to MercuryDPM.  Comput. Phys. Commun. 249, 107129 (2020).  DOI 10.1016/j.cpc.2019.107129
  * (22) Lubachevsky, B. D., Stillinger, F. H.: Geometric properties of random disk packings.  J. Stat. Phys. 60, 561–583 (1990).  DOI 10.1007/BF01025983