1 Introduction

In this paper, we consider a storage system with Lévy netput. In other words, the workload process \(\{Q(t), t\ge 0\}\) is a spectrally one-sided Lévy process X(t) that is reflected at 0:

$$\begin{aligned} Q(t):=X(t)-{{\text {inf}}_{s\le t} }(X(s))^-, \end{aligned}$$
(1)

where \({A^{-}}=\min \{A,0\}\). We assume that the drift of the process X(t) is negative; that is, we have \({\mathbb E}X(1)<0\). This stability condition guarantees the existence of a stationary distribution \(\pi \) of Q, which by virtue of “Reich’s identity” can be expressed in terms of the all-time supremum:

$$\begin{aligned} \pi (x)={\mathbb P}\left( \sup _{t\ge 0} X(t)\le x\right) . \end{aligned}$$
(2)

In the sequel, we consider the initial distribution Q(0) sampled from this steady-state distribution, which is indicated by adding the subscript \(\pi \) to the probability measure \({\mathbb P}\) and to the associated expectation \({\mathbb E}\).

Now, let T denote the busy period; that is,

$$\begin{aligned} T=\inf \{t\ge 0:\ Q(t)=0\}.\end{aligned}$$

We will further consider the Yaglom limit

$$\begin{aligned} \mu (\mathrm{d} x, \mathrm{d} y):= \lim _{t\rightarrow \infty }{\mathbb P}_\pi (Q(0)\in d x, Q(t)\in \mathrm{d} y\,|\,T >t), \end{aligned}$$

where the convergence is to be understood in the weak sense. Yaglom [30] limits form a probability measure and a particular case of the quasi-stationary (QS) distribution, which is an invariant distribution for the process conditioned on non-extinction; that is, we condition on the event that the process survives some killing event (for example, related to exiting from some subset of possible values).

Yaglom [30] was the first to explicitly identify QS distributions for the subcritical Bienaymé–Galton–Watson branching process. This result has been generalized in the context of the continuous-time branching process and the Fleming–Viot process; see [1, 9, 18]. Similar results were also derived for Markov chains on positive integers with an absorbing state at the origin; see Seneta and Vere-Jones [27], Tweedie [29], Jacka and Roberts [13] and the bibliographic database of Pollet [25]. Recently, Foley and McDonald showed that Yaglom limits may depend on the starting state [10].

The research on QS distributions has been very extensive. Martinez and San Martin [22] analyze the Brownian motion with drift exiting from the positive half-line, complementing the result for random walks obtained by Iglehart [12]. In addition, QS laws have been studied for various Lévy processes. Kyprianou [15] found the Laplace transform of the QS distribution for the workload process of the stable M/G/1 queue with service times that have a rational moment generating function. Kyprianou and Palmowski [17] identified the QS distribution associated with a general light-tailed Lévy process. Haas and Rivero [26] found (after appropriate scaling) the QS distribution when the Lévy process under study has a jump measure with a regularly varying tail. The speed of convergence (in total variation) to the QS distribution for population processes has been studied in [5]. Finally, Mandjes et al. [21] derived the QS distribution of the workload process Q(t). This paper builds upon [21].

A contribution of this paper lies in proving that the speed of convergence to the quasi-stationary distribution is surprisingly slow (of order 1/t). We also identify a measure \(\xi (dx,dy)\) (which we call the second-order quasi-stationary measure) such that

$$\begin{aligned} \lim _{t\rightarrow \infty }t\times |{\mathbb P}_\pi (Q(0)\in \mathrm{d} x, Q(t)\in \mathrm{d} y\,|\,T >t)-\mu (\mathrm{d} x, \mathrm{d} y)|=\xi (\mathrm{d}x,\mathrm{d}y), \end{aligned}$$
(3)

where the above convergence is in the weak sense. We hence prove the conjecture posed in Polak and Rolski [24], which proved the above statement for a birth–death process by using an asymptotic expansion of a transition function and certain properties of Bessel functions. In this paper, we suggest a new method, which relies on a refined Tauberian-type expansion of the Laplace transform. We also analyze in detail the \(M/E(2,\nu )/1\) queue and a Brownian-driven queue.

If we want to simulate the quasi-stationary distribution directly from its definition given in (3), then the result stated in (3) shows that the speed of such a simulation is very slow. Still, in the literature there are papers giving other efficient algorithms of simulation of quasi-stationary measures; see, for example, Blanchet et al. [4] and references therein.

The main result in (3) contrasts the typical results derived for the regular stationary distribution of Markov processes where, in most of the cases, the rate is exponential. More precisely, for many models the distance between the distribution of the stochastic process at time t and its stationary distribution decays exponentially fast in t. The typical distances used are the total variation distance, the separation distance, and the \(L^2\) distance. The classical results concern mainly Markov chains and use Perron and Frobenius theory, renewal equations or the coupling method; see, for example, [6, 8, 14, 19, 20] and references therein. Another method concerns Harris recurrent Markov processes, and it is based on the construction of a special Lyapunov function and then the application of Foster-Lyapunov criteria; see, for example, [2, 23, 28]. All the above-mentioned methods, though, are different from the one used in this paper, which is based on expansions of Laplace transforms.

The paper is organized as follows: In next section, we introduce the notation and basic facts that are used later. In Sect. 3, we present the main results. The central step for the proof of the main results is given in Sect. 4. Finally, the last section provides some examples.

2 Preliminaries

We follow [16] for definitions, notation and basic facts on Lévy processes. Let \(X\equiv (X(t))_t\) be a Lévy process, which is defined on the filtered space \((\Omega ,\mathcal { F},\{\mathcal {F}_t\}_{t\ge 0}, {\mathbb P})\) with the natural filtration that satisfies the usual assumptions of right continuity and completion. We define \({\mathbb P}_x\) as \({\mathbb P}_x(X(0)=x)=1\) and \({\mathbb P}_0={\mathbb P}\); similarly, \({\mathbb E}_x\) is the expectation with respect to \({\mathbb P}_x\). We denote by \(\Pi (\cdot )\) the jump measure of X. If X is spectrally negative (resp. spectrally positive), then \(\Pi \) is supported in the non-positive half-line (resp. in the non-negative half-line); in other words, jumps are non-positive (resp. non-negative). We define the Laplace exponent \(\psi ({\eta })\) by

$$\begin{aligned} {\mathbb E}e^{{\eta } X(t)}=e^{t\psi ({\eta })}, \end{aligned}$$
(4)

for \({\eta }\in \mathbb {R}\) such that the left-hand side of (4) is well defined (which holds at least for \(\eta \ge 0\)). We denote by \(\Phi (s):=\inf \{{\eta }\ge 0: \psi ({\eta })>s\}\) the right inverse of \(\psi \); see [16] for details.

Dual process We also consider the dual process \(\hat{X}_t=-X_t\) with jump measure \(\hat{\Pi }\left( 0,y\right) =\Pi \left( -y,0\right) \). The characteristics of \(\hat{X}\) are indicated by using the same symbols as for X, but with a ‘\(\hat{}\)’ added. We will write

$$\begin{aligned} \hat{\psi }(\eta )= t^{-1} \log {\mathbb E}e^{{\eta } \hat{X}(t)}=t^{-1} \log {\mathbb E}e^{-{\eta } X(t)}=\psi (-\eta ). \end{aligned}$$
(5)

We skip the symbol ‘\(\hat{}\)’ for Q(t) and hence for T and all quantities related to Q, as it will be clear from the context if a statement concerns the spectrally negative or the spectrally positive case.

Asymptotic expansions

Consider a function \(f:\mathbb {R}\rightarrow \mathbb {R}\) such that \(f(x)=0\) for \(x<0\). Let \(\tilde{f}(z):=\int _0^\infty e^{-z x}f(x)\,d x\) for \(z\in \mathbb {C}\) be its Laplace transform. Consider singularities of \(\tilde{f}(z)\); among these, let \(\mathfrak {R}(a_0)<0\) be the one with the largest real part. Notice that this yields the integrability of \(\int _0^\infty |f(x)|\, d x\). The inversion formula reads

$$\begin{aligned} f(x)=\frac{1}{2\pi \mathrm{i}}\int _{a-\mathrm{i}\infty }^{a+\mathrm{i}\infty }\tilde{f}(z) e^{zx}\, d z \end{aligned}$$

for some (and then any) \(\mathfrak {R}(a)>\mathfrak {R}(a_0)\). In this paper, we need the following Tauberian theorem, found in Doetsch [7, Theorem 37.1], where the behaviour of the Laplace transform around the singularity \(a_0\) plays a crucial role.

First, recall the concept of the \(\mathfrak W\)-contour, centered at \(a_0\), with a half-angle of opening \(\pi /2<\psi \le \pi \), as depicted in [7, Fig. 30, p. 240]; also, for the purposes of our problem, \({\mathcal G}_{a_0}(\psi )\) is the region between the contour \(\mathfrak W\) and the line \(\mathfrak {R}(z)=0\). More precisely,

$$\begin{aligned} {\mathcal G}_{a_0}(\psi ) \equiv \{z \in \mathbb {C}; \mathfrak {R}(z)< 0, z \ne a_0, |\arg (z - a_0)| < \psi \}, \end{aligned}$$

where \(\arg z\) is the principal part of the argument of the complex number z. In the following theorem, conditions are identified that provide an asymptotic expansion of the Laplace transform.

Theorem 1

([7, Theorem 37.1]) Suppose that for \(\tilde{f}:\mathbb {C}\rightarrow \mathbb {C}\) and \(\mathfrak {R}(a_0) < 0\) the following three conditions hold:

(A1):

\(\tilde{f}(\cdot )\) is analytic in a region \({\mathcal G}_{a_0}(\psi )\) for some \(\pi /2<\psi \le \pi \);

(A2):

\(\tilde{f}(z) \rightarrow 0\) as \(|z| \rightarrow \infty \) for \(z \in {\mathcal G}_{a_0}(\psi )\);

(A3):

for some constants \(c_\nu \), \(\tilde{f}(s)\) has in \(|\mathrm {arg} (s - a_0)| < \psi \) the asymptotic expansion

$$\begin{aligned} \tilde{f}(s)\approx \sum _{\nu =0}^\infty c_\nu (s-a_0)^{\lambda _\nu }, \qquad (\mathfrak {R}(\lambda _0)<\mathfrak {R}(\lambda _1)<\ldots )\text {\ as\ }s\rightarrow a_0. \end{aligned}$$
(6)

Then we conclude that, as \(t\rightarrow \infty \), f(t) has the asymptotic expansion

$$\begin{aligned} f(t)\approx e^{a_0 t}\sum _{\nu =0}^\infty \frac{c_\nu }{\Gamma (-\lambda _\nu )}\frac{1}{t^{\lambda _\nu +1}},\qquad \left( \frac{1}{\Gamma (-\lambda _\nu )}=0\text {\ for\ } \lambda _\nu =0,1,2,\ldots \right) . \end{aligned}$$

Assumptions For a spectrally negative Lévy process X(t), we impose the following assumptions:

(SN)

  1. (SN1)

    \(\psi ({\vartheta })\) attains its strictly negative minimum at \({\vartheta }{^*}>0\) (and hence \(\psi '({\vartheta }{^*})=0\));

  2. (SN1)

    \(\Phi \) is analytic in \(\mathcal {G}_{\zeta {^*}}(\phi )\) for \( \pi /2< \phi \le \pi \), where

    $$\begin{aligned} \zeta {^*}:=\psi ({\vartheta }{^*})<0. \end{aligned}$$
    (7)

Similar conditions are assumed for a spectrally positive Lévy process X(t), for which \(\hat{X}_t\) is spectrally negative with Laplace exponent \(\hat{\psi }\) :

(SP)

  1. (SP1)

    \(\hat{\psi }({\vartheta })\) attains its strictly negative minimum at \({\vartheta }{^*}<0\) (and hence \(\hat{\psi }'({\vartheta }{^*})=0\));

  2. (SP1)

    \(\hat{\Phi }\) is analytic in \(\mathcal {G}_{\zeta {^*}}(\phi )\) for \( \pi /2< \phi \le \pi \), where

    $$\begin{aligned} \zeta {^*}:=\hat{\psi }({\vartheta }{^*})<0. \end{aligned}$$
    (8)

To check the above assumptions, we can use the concept of semiexponentiality of a function f (see [11, p. 314]).

Definition 2

(Semiexponentiality) A function f is said to be semiexponential if, for some \(0< \phi \le \pi /2\), there exists a finite and strictly negative function \(\gamma ({\vartheta })\), called the indicator function, defined as the infimum of all \(a\in \mathbb {R}\) such that

$$\begin{aligned}\left| f(e^{\mathrm{i}{\vartheta }}r)\right| <e^{ar}\end{aligned}$$

for all sufficiently large r; here \(-\phi \le {\vartheta }\le \phi \) and \(\sup \gamma ({\vartheta })<0\).

It was proved in [21] using [11, Thm. 10.9f] that if there exists a density of \(\Pi \) (resp. \(\hat{\Pi }\)) which is of semiexponential type, then \(\Phi \) (resp. \(\hat{\Phi }\)) is analytic in \(\mathcal {G}_{\zeta {^*}}(\phi )\) for \( \pi /2< \phi \le \pi \). In particular, this assumption holds, for example, for a Brownian motion \(X(t)=\sigma B(t)-c t\) with linear drift, where \(c>0\) and B is the standard Brownian motion.

Quasi-stationary distribution Denote by

$$\begin{aligned}\tilde{\mu }(\alpha , \beta ):=\int _0^\infty \int _0^\infty e^{-\alpha x}e^{-\beta y}\mu (dx,dy)\end{aligned}$$

the bivariate Laplace transform of the quasi-stationary measure. The quasi-stationary distribution of the workload process Q(t) in the stationary regime was identified in [21].

Theorem 3

(Mandjes et al. [21])

(i) Under (SN),

$$\begin{aligned} \tilde{\mu }(\alpha ,\beta )= \frac{-\psi ({\vartheta }{^*})}{\psi (\alpha +\Phi (0))-\zeta {^*}}\frac{({\vartheta }{^*})^2}{({\vartheta }{^*}+\beta )^2}\;. \end{aligned}$$

(ii) Under (SP),

$$\begin{aligned} \tilde{\mu }(\alpha ,\beta )= \hat{\psi }^2({\vartheta }{^*})\cdot \frac{\hat{\psi }(\alpha +{\vartheta }{^*})-(\alpha +{\vartheta }{^*})\hat{\psi }'(\alpha +{\vartheta }{^*})}{\hat{\psi }^2(\alpha +{\vartheta }{^*})(\hat{\psi }({\vartheta }{^*}) -\hat{\psi }(\beta ))}\;. \end{aligned}$$

Note that conditions (SN) and (SP) are satisfied since we assumed that \({\mathbb E}X_1<0\). The two key components of the proof of this result are based first on Wiener–Hopf factorization, from which master formulas can be derived (given below), and second on either some expansion theorems (see [3] and [17]) or on some Tauberian-type theorems.

Master formulas Recall that Q(t) given in (1) is a workload process with stationary distribution (2) and busy period T. We define now the double Laplace–Stieltjes transform:

$$\begin{aligned} L({\vartheta };\alpha ,\beta ):=\int _0^\infty e^{-{\vartheta } t} {\mathbb E}_\pi [e^{-\alpha Q(0)-\beta Q(t)},T>t] \, d t. \end{aligned}$$

In [21], the following representations of L were derived.

Proposition 4

[Mandjes et al. [21]]

  1. (i)

    Under (SN),

    $$\begin{aligned} L({\vartheta };\alpha ,\beta )= \frac{\Phi ({\vartheta })-\alpha -\Phi (0)}{\Phi ({\vartheta })+\beta }\frac{\Phi (0)}{\alpha +\beta +\Phi (0)}\frac{1}{{\vartheta }-\psi (\alpha +\Phi (0))}. \end{aligned}$$
  2. (ii)

    Under (SP),

    $$\begin{aligned} L({\vartheta };\alpha ,\beta )=\frac{\hat{\psi }'(0+)}{{\vartheta }-\hat{\psi }(\beta )}\left( \frac{\alpha +\beta }{\hat{\psi }(\alpha +\beta )} - \frac{\alpha +\hat{\Phi }({\vartheta })}{\hat{\psi }(\alpha +\hat{\Phi }({\vartheta }))}\right) . \end{aligned}$$

From the proposition above, it follows that under assumptions (SN) or (SP) one can extend analytically \(L({\vartheta };\alpha ,\beta )\) into \({\mathcal G}_{\zeta {^*}}(\psi )\) for some \(\pi /2<\psi \le \pi \).

To prove the main result we will expand \(L(\vartheta , \alpha ,\beta )\) around \(\vartheta =\zeta ^*\) using Theorem 3. Then we will apply Theorem 1 to identify the expansion of \({\mathbb E}_\pi [e^{-\alpha Q(0)-\beta Q(t)}|T>t]\) for large values of t. In the last step we will use the definition of the quasi-stationary distribution.

3 Main results

We state now the main results of this paper. Define the constants

$$\begin{aligned} A_1&:=\sqrt{\frac{2}{ \psi {''}(\vartheta ^*)}},&B_1&:=\sqrt{\frac{2}{ \hat{\psi }{''}(\vartheta ^*)}},\\ A_2&:=-\frac{\psi {'''}(\vartheta ^*)}{3\psi {''}(\vartheta ^*)^2},&B_2&:=-\frac{\hat{\psi }{'''}(\vartheta ^*)}{3\hat{\psi }{''}(\vartheta ^*)^2},\\ A_3&:=\frac{5 \psi ^{(3)}\left( \theta ^*\right) ^2-3 \psi ^{(4)}\left( \theta ^*\right) \psi ''\left( \theta ^*\right) }{18 \sqrt{2} \psi ''\left( \theta ^*\right) ^{7/2}}&B_3&:=\frac{5 \hat{\psi }^{(3)}\left( \theta ^*\right) ^2-3 \hat{\psi }^{(4)}\left( \theta ^*\right) \hat{\psi }''\left( \theta ^*\right) }{18 \sqrt{2} \hat{\psi }''\left( \theta ^*\right) ^{7/2}}. \end{aligned}$$

We start with the following expansion for the double Laplace–Stieltjes transform.

Proposition 5

If (SN) or (SP) hold, then

$$\begin{aligned} L({\vartheta };\alpha ,\beta ) {=}C_0(\alpha ,\beta ){+} C_1(\alpha ,\beta ) (\vartheta {-}\zeta ^*)^{1/2} {+} C_2(\alpha ,\beta ) ({\vartheta }{-}\zeta ^{*}) {+} C_3(\alpha ,\beta ) (\vartheta {-}\zeta ^{*})^{3/2} {+} o((\vartheta {-}\zeta ^{*})^{3/2}), \end{aligned}$$

for \(\zeta ^*<0\) defined in (7) and (8).

(i):

Under (SN),

$$\begin{aligned} C_0(\alpha ,\beta )&= -\frac{\Phi (0) }{\alpha +\beta +\Phi (0)} \frac{\alpha -\vartheta ^*+\Phi (0)}{\left( \beta +\vartheta ^*\right) \left( \zeta ^*-\psi (\alpha +\Phi (0))\right) },\\ C_1(\alpha ,\beta )&= \frac{A_1 \Phi (0)}{\left( \beta +\vartheta ^*\right) ^2 \left( \zeta ^*-\psi (\alpha +\Phi (0))\right) },\\ C_2(\alpha ,\beta )&= \frac{\Phi (0)\left( \frac{\left( \beta +\vartheta ^*\right) ^2 \left( \alpha -\vartheta ^*+\Phi (0)\right) }{\alpha +\beta +\Phi (0)}-\left( A_2 \left( \beta +\vartheta ^*\right) -A_1^2\right) \left( \psi (\alpha +\Phi (0))-\zeta ^*\right) \right) }{\left( \beta +\vartheta ^*\right) ^3\left( \zeta ^*-\psi (\alpha +\Phi (0))\right) ^2},\\ C_3(\alpha ,\beta )&=\frac{\Phi (0) A_3}{\left( \beta +\vartheta ^*\right) ^2\left( \zeta ^*-\psi (\alpha +\Phi (0))\right) } -\frac{\Phi (0) }{\left( \beta +\vartheta ^*\right) ^4\left( \zeta ^*-\psi (\alpha +\Phi (0))\right) ^2}\\&\quad \cdot \left( A_1 \left( \beta +\vartheta ^*\right) \left( 2 A_2 \left( \zeta ^*-\psi (\alpha +\Phi (0))\right) +\beta +\vartheta ^*\right) \right. \\&\qquad \qquad \left. +A_1^3\left( \psi (\alpha +\Phi (0))-\zeta ^*\right) \right) . \end{aligned}$$
(ii):

Under (SP),

$$\begin{aligned} C_0(\alpha ,\beta )&= \left( \frac{\alpha +\beta }{\hat{\psi }(\alpha +\beta )}-\frac{\alpha +\vartheta ^*}{\hat{\psi }\left( \alpha +\vartheta ^*\right) }\right) \frac{\hat{\psi }'(0)}{\zeta ^*-\hat{\psi }\left( \beta \right) }, \\ C_1(\alpha ,\beta )&= -\hat{\psi }'(0)B_1\frac{ \hat{\psi }\left( \alpha +\vartheta ^*\right) -(\alpha +\vartheta ^*) \hat{\psi }'\left( \alpha +\vartheta ^*\right) }{\hat{\psi }\left( \alpha +\vartheta ^*\right) ^2 \left( \zeta ^*-\hat{\psi }(\beta )\right) }, \\ C_2(\alpha ,\beta )&=\frac{ \hat{\psi }'(0)}{\zeta ^*-\hat{\psi }(\beta )} \Bigg [ \frac{B_1^2 \hat{\psi }'\left( \alpha +\vartheta ^*\right) }{\hat{\psi }\left( \alpha +\vartheta ^*\right) ^2} -\frac{B_2}{\hat{\psi }\left( \alpha +\vartheta ^*\right) } -\frac{\frac{\alpha +\beta }{\hat{\psi }(\alpha +\beta )}-\frac{\alpha +\vartheta ^*}{\hat{\psi }\left( \alpha +\vartheta ^*\right) }}{\zeta ^*-\hat{\psi }(\beta )} \\&\quad +\frac{\left( \alpha +\vartheta ^*\right) \left( B_1^2 \left( \hat{\psi }\left( \alpha +\vartheta ^*\right) \hat{\psi }''\left( \alpha +\vartheta ^*\right) -2 \hat{\psi }'\left( \alpha +\vartheta ^*\right) ^2\right) +2 B_2 \hat{\psi }\left( \alpha +\vartheta ^*\right) \hat{\psi }'\left( \alpha +\vartheta ^*\right) \right) }{2 \hat{\psi }\left( \alpha +\vartheta ^*\right) ^3} \Bigg ], \end{aligned}$$
$$\begin{aligned} C_3(\alpha ,\beta )&=\hat{\psi }'(0) \Bigg [B_1\frac{\hat{\psi }(\alpha +\vartheta ^*)-\hat{\psi }'(\alpha +\vartheta ^*)(\alpha + \vartheta ^*)}{\hat{\psi }(\alpha +\vartheta ^*)^2 \left( \zeta ^*-\hat{\psi }(\beta )\right) ^2} \\&\quad -\frac{1}{6 \hat{\psi }\left( \alpha +\vartheta ^*\right) ^4 \left( \zeta ^*-\hat{\psi }(\beta )\right) }\Bigg (6 B_3 \hat{\psi }\left( \alpha +\vartheta ^*\right) ^3\\&\quad -6 {B_{1}}^{3} \left( \alpha +\vartheta ^*\right) \hat{\psi }'\left( \alpha +\vartheta ^*\right) ^3 \\&\quad +6 B_1 \hat{\psi }\left( \alpha +\vartheta ^*\right) \hat{\psi }'\left( \alpha +\vartheta ^*\right) \\&\qquad \left( B_1^2 \left( \alpha +\vartheta ^*\right) \hat{\psi }''\left( \alpha +\vartheta ^*\right) +\hat{\psi }'\left( \alpha +\vartheta ^*\right) \left( B_1^2+2 B_2 (\alpha +\vartheta ^*)\right) \right) \\&\quad -6 \hat{\psi }(\alpha +\vartheta ^*)^2\hat{\psi }'\left( \alpha +\vartheta ^*\right) \left( 2 B_1 B_2+ B_3 (\alpha +\vartheta ^*)\right) \\&\quad -\hat{\psi }(\alpha +\vartheta ^*)^2 B_1 \left( 3 \hat{\psi }''\left( \alpha +\vartheta ^*\right) \left( B_1^2+2 B_2(\alpha +\vartheta ^*)\right) \right. \\&\quad \left. + B_1^2 \left( \alpha +\vartheta ^*\right) \hat{\psi }^{(3)}\left( \alpha +\vartheta ^*\right) \right) \Bigg )\Bigg ]. \end{aligned}$$

Theorem 6

If assumptions (SN) hold for a spectrally negative Lévy process X or (SP) hold for a spectrally positive Lévy process \(\hat{X}\), then the measure \(\xi (dx,dy)\) defined formally in (3) exists. That is, the speed of convergence to the quasi-stationary distribution is of order 1/t. Moreover, the Laplace transform \(\tilde{\xi } (\alpha , \beta ):=\int _0^\infty \int _0^\infty e^{-\alpha x}e^{-\beta y}\xi (dx,dy)\) equals

$$\begin{aligned} \tilde{\xi } (\alpha , \beta )=\frac{(-3)}{2C_1(0,0)}\left( C_3(\alpha ,\beta )-\tilde{\mu }(\alpha ,\beta )C_3(0,0)\right) . \end{aligned}$$
(9)

Proof

Recall that under the imposed assumptions (SN) or (SP), the Laplace exponent \(L({\vartheta };\alpha ,\beta )\) as a function of \(\vartheta \) satisfies the assumptions of Theorem 1. Hence, from Proposition 5 and Theorem 1 we have that, as \(t\rightarrow \infty \),

$$\begin{aligned}&{\mathbb E}_\pi [e^{-\alpha Q(0)-\beta Q(t)},T>t]\nonumber \\&\quad =e^{\zeta ^*t}\left( \frac{C_1(\alpha ,\beta )}{\Gamma (-1/2)}t^{-3/2}+\frac{C_3(\alpha ,\beta )}{\Gamma (-3/2)}t^{-5/2}+o(t^{-5/2})\right) . \end{aligned}$$
(10)

Further, note that by Theorem 3

$$\begin{aligned} \tilde{\mu }(\alpha ,\beta )=C_1(\alpha ,\beta )/C_1(0,0). \end{aligned}$$
(11)

Now, straightforward calculations give

$$\begin{aligned}&{\mathbb E}_\pi [e^{-\alpha Q(0)-\beta Q(t)}|T>t]\\&\quad =\frac{C_1(\alpha ,\beta )}{C_1(0,0)} +\frac{\Gamma (-1/2)}{\Gamma (-3/2)}\left( \frac{C_3(\alpha ,\beta )}{C_1(0,0)}-\frac{C_1(\alpha ,\beta )C_3(0,0)}{C_1^2(0,0)}\right) t^{-1} + o(t^{-1}), \end{aligned}$$

which completes the proof due to the definition of the measure \(\xi \).

Remark 7

Formally starting from Proposition 4, then using a generalization of Proposition 5 in the next step, combined with Theorem 1, will give more terms of the expansion of \({\mathbb E}_\pi [e^{-\alpha Q(0)-\beta Q(t)}|T>t]\), and hence a longer expansion of the measure \({\mathbb P}_\pi (Q(0)\in dx, Q(t)\in dy|T>t)\) as \(t\rightarrow \infty \).

4 Proof of Proposition 5

Presume now that assumptions (SP) hold. We start from deriving the expansion of \(\hat{\Phi }(\vartheta )\):

$$\begin{aligned} \hat{\Phi }(s)=\vartheta ^* + B_1 (s-\zeta ^{*})^{1/2} + B_2 ({s}-\zeta ^{*}) + B_3 (s-\zeta ^{*})^{3/2} + o((s-\zeta ^{*})^{3/2}) \end{aligned}$$
(12)

as \(s\downarrow \zeta ^*\). Indeed, from a Taylor series expansion of \(\hat{\psi }\) around \(\vartheta ^*\) and the condition that \(\hat{\psi }'({\vartheta ^*})=0\), we have

$$\begin{aligned}&\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})\nonumber \\&\quad =\frac{({\vartheta }-{\vartheta ^*})^2}{2} \hat{\psi }''({\vartheta ^*})\nonumber \\&\qquad +\frac{({\vartheta }-{\vartheta ^*})^3}{6} \hat{\psi }^{(3)}({\vartheta ^*})\nonumber \\&\qquad +\frac{({\vartheta }-{\vartheta ^*})^4}{24} \hat{\psi }^{(4)}({\vartheta ^*})+ o(({\vartheta }-{\vartheta ^*})^4). \end{aligned}$$
(13)

From the above equation we can derive

$$\begin{aligned} \vartheta -\vartheta ^*= & {} B_1 \sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}+ B_2\left( \hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})\right) +B_3 \left( \hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})\right) ^{3/2}\nonumber \\&\quad +\,o(({\vartheta }-\vartheta ^*)^{3/2}). \end{aligned}$$
(14)

The main idea of getting (14) is to include this expansion (14) into (13) and to match respective powers of \(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})\). Indeed, let us start from the first-order expansion

$$\begin{aligned} \hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})= \frac{({\vartheta }-{\vartheta ^*})^2}{2} \hat{\psi }''({\vartheta ^*}) +o(({\vartheta }-{\vartheta ^*})^2)=\frac{({\vartheta }-{\vartheta ^*})^2}{2} \hat{\psi }''({\vartheta ^*})(1+o(1)).\nonumber \\ \end{aligned}$$
(15)

Then

$$\begin{aligned} \vartheta -\vartheta ^*&=B_1\sqrt{(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*}))/(1+o(1))}= B_1\sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}\sqrt{1+o(1)} \nonumber \\&=B_1\sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}(1+o(1))\nonumber \\&=B_1\sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}+o\left( \sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}\right) \end{aligned}$$
(16)
$$\begin{aligned}&=B_1\sqrt{(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*}))} +o(\vartheta -\vartheta ^*). \end{aligned}$$
(17)

Note that in the last equation \(o\left( \sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}\right) = o(\vartheta -\vartheta ^*)\) by (15). Now, we continue with the second-order computation:

$$\begin{aligned}\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})= \frac{({\vartheta }-{\vartheta ^*})^2}{2} \hat{\psi }''({\vartheta ^*})+ \frac{({\vartheta }-{\vartheta ^*})^3}{6} \hat{\psi }^{(3)}({\vartheta ^*})+o\left( ({\vartheta }-{\vartheta ^*})^3\right) .\end{aligned}$$

Then putting expression (16) in all powers of \({\vartheta }-{\vartheta ^*}\) in the increments on the right-hand side of the above equation, we can conclude that the last increment in (16) equals

$$\begin{aligned}o\left( \sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}\right)&= \frac{-\left( B_1\sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}\right) ^3}{6 B_1\sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}\hat{\psi }''({\vartheta ^*})} \hat{\psi }^{(3)}({\vartheta ^*})\\&\quad +o\left( ({\vartheta }-{\vartheta ^*})^3\right) +o(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*}))\\&=B_2(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})) +o(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})), \end{aligned}$$

where the last equality holds by observing that \(o\left( ({\vartheta }-{\vartheta ^*})^3\right) = o(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*}))\), which follows from (15). Hence, plugging this back into (16) we derive

$$\begin{aligned} \vartheta -\vartheta ^*=B_1\sqrt{\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})}+B_2(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})) +o(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*})). \end{aligned}$$
(18)

In the last step, we do the third-order computations required to identify \(o(\hat{\psi }({\vartheta })-\hat{\psi }({\vartheta ^*}))\) in (18) by again plugging (18) to (13). This produces (14). Substituting \({\vartheta }=\hat{\Phi }(s)\) and using \(\hat{\psi }(\hat{\Phi }(s))=s\) completes the proof of (12).

Having proven (12), we plug it into Proposition 4 and order the outcome according to powers of \(s-\zeta ^{*}\). This will complete the proof of the spectrally positive case (SP). Similarly, when (SN) holds then

$$\begin{aligned} {\Phi }(s)=\vartheta ^* + A_1 (s-\zeta ^*)^{1/2} + A_2 ({s}-\zeta ^{*}) + A_3 (s-\zeta ^{*})^{3/2} + o((s-\zeta ^{*})^{3/2}) \end{aligned}$$

as \(s\downarrow \zeta ^*\). Then, using Proposition 4 in the same way as before gives the required assertion after some simple manipulations. \(\square \)

5 Examples

In this section, we illustrate our theory through two examples.

Example 8

In this case

$$\begin{aligned} X(t) = \sum _{i=1}^{N(t)} \sigma _i-t, \end{aligned}$$
(19)

where \(\sigma _i\) (where \(i=1, 2, \ldots \)) are i.i.d. service times that have an Erlang\((2,\nu )\) distribution, i.e., with rate \(\nu \) and two phases. The arrival process is a homogeneous Poisson process N(t) with rate \(\lambda \). We assume that \(\varrho :=2\lambda /\nu <1\). For the Laplace exponent, we have that

$$\begin{aligned}\hat{\psi }(\eta )=\eta -\lambda +\lambda \left( \frac{\nu }{\eta + \nu }\right) ^2,\end{aligned}$$

which attains its minimum at

$$\begin{aligned} {\vartheta }^*=\root 3 \of {2\lambda \nu ^2} -\nu \end{aligned}$$

and it is equal to

$$\begin{aligned} \zeta ^*=\frac{3\root 3 \of {\lambda \nu ^2}}{\root 3 \of {4}}-\nu -\lambda .\end{aligned}$$

One can easily check that all assumptions (SP) are satisfied. In particular, \(\hat{\Phi }(z)\) is analytic in \(\mathbb {C}\setminus (-\infty ,\zeta {^*}]\). Then, Proposition 5 gives for \(C_1(\alpha ,\beta )\) that

$$\begin{aligned}&C_1(\alpha ,\beta )\\&\quad =\frac{2^{5/3} \lambda (\beta +\nu )^2 (2 \lambda -\nu ) \sqrt{\root 3 \of {\lambda } \nu ^{2/3}} \left( \alpha +\root 3 \of {2} \root 3 \of {\lambda } \nu ^{2/3}\right) \left( \alpha +\root 3 \of {2} \root 3 \of {\lambda } \nu ^{2/3}+2 \nu \right) }{\sqrt{3} \nu \left( \alpha ^2+2 \root 3 \of {2} \alpha \root 3 \of {\lambda } \nu ^{2/3}-\lambda (\alpha +\nu )+2^{2/3} \lambda ^{2/3} \nu ^{4/3}-\root 3 \of {2} \lambda ^{4/3} \nu ^{2/3}\right) ^2}\ \\&\qquad \times \, \frac{1}{\left( 2 \beta ^3-3 \root 3 \of {2} \beta ^2 \root 3 \of {\lambda } \nu ^{2/3}+6 \beta ^2 \nu -6 \root 3 \of {2} \beta \root 3 \of {\lambda } \nu ^{5/3}+2 \nu ^2 (3 \beta +\lambda )-3 \root 3 \of {2} \root 3 \of {\lambda } \nu ^{8/3}+2 \nu ^3\right) }, \end{aligned}$$

which is sufficient to compute the bivariate Laplace transform of the quasi-stationary measure given in (11). As a numerical illustration, if \(\lambda =1\) and \(\nu =4\), then

$$\begin{aligned}&\tilde{\mu }(\alpha ,\beta )\\&\quad =\frac{C_1(\alpha ,\beta )}{C_1(0,0)}=\bigg \{\left( 10-3*2^{5/3}\right) ^2 \left( \alpha +2^{5/3}\right) \left( \alpha +2^{5/3}+8\right) (\beta +4)^2\bigg \}\\&\qquad \bigg / \bigg \{2 \left( 2 \beta ^3-6*2^{2/3} \beta ^2+24 \beta ^2-3*2^{14/3} \beta +32(3 \beta +1)-3*2^{17/3}+2^7\right) \\&\qquad \times \, \left( \alpha ^2+(2^{8/3}-1)\alpha -2^{5/3}+2^{10/3}-4\right) ^2\bigg \}. \end{aligned}$$

To identify the bivariate Laplace transform of the second-order quasi-stationary measure \(\xi \) given in Theorem 6, we need to find \(C_3(\alpha , \beta )\). Its expression is rather complex:

$$\begin{aligned}&C_3(\alpha ,\beta )=\frac{1-\frac{2 \lambda }{\nu }}{6 \left( -\frac{\lambda \nu ^2}{(\beta +\nu )^2}-\beta +\frac{3 \root 3 \of {\lambda } \nu ^{2/3}}{2^{2/3}}-\nu \right) \left( \frac{\lambda \nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^2}+\alpha +\theta ^*-\lambda \right) ^4}\\&\quad \times \Bigg \{\frac{2^{5/3} \sqrt{3} \lambda \left( \alpha +\theta ^*\right) ^2 \sqrt{\root 3 \of {\lambda } \nu ^{2/3}} \left( \alpha +\theta ^*+3 \nu \right) \left( \lambda \left( \frac{\nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^2}-1\right) +\alpha +\theta ^*\right) ^2}{\left( \alpha +\theta ^*+\nu \right) ^3 \left( \frac{\lambda \nu ^2}{(\beta +\nu )^2}+\beta -\frac{3 \root 3 \of {\lambda } \nu ^{2/3}}{2^{2/3}}+\nu \right) }\\&\quad +\frac{2}{3^{7/3}}\Bigg [36 \left( \alpha +\theta ^*\right) \left( \root 3 \of {\lambda } \nu ^{2/3}\right) ^{3/2} \left( 1-\frac{2 \lambda \nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^3}\right) ^3 \\&\quad -27* 2^{2/3} \sqrt{\root 3 \of {\lambda } \nu ^{2/3}}\left( 1-\frac{2 \lambda \nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^3}\right) \\&\quad \times \bigg [\frac{{2}^{7/2} \lambda ^{4/3} \nu ^{8/3} \left( \alpha +\theta ^*\right) }{\left( \alpha +\theta ^*+\nu \right) ^4}+\Big [ \frac{8\left( \alpha +\theta ^*\right) }{9} +\frac{2}{3} \root 3 \of {2} \root 3 \of {\lambda } \nu ^{2/3} \left( 1-\frac{2 \lambda \nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^3}\right) \Big ] \bigg ]\\&\quad \times \, \left( \lambda \left( \frac{\nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^2}-1\right) +\alpha +\theta ^*\right) +\root 3 \of {\frac{2}{\lambda \nu ^2}} \left( 5 (\alpha +\theta ^*)+3 \root 3 \of {2^{10}\lambda \nu ^{2}}\right) \left( 1-\frac{2 \lambda \nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^3}\right) \\&\quad \times \left( \lambda \left( \frac{\nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^2}-1\right) +\alpha +\theta ^*\right) ^2-\frac{5 \root 3 \of {2} \left( \lambda \left( \frac{\nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^2}-1\right) +\alpha +\theta ^*\right) ^3}{\lambda \nu ^2 \left( \frac{1}{\root 3 \of {\lambda } \nu ^{2/3}}\right) ^{5/2}} + \frac{9\ 2^{2/3} \lambda \nu ^2 \sqrt{\root 3 \of {\lambda } \nu ^{2/3}}}{\left( \alpha +\theta ^*+\nu \right) ^5}\\&\quad \times \, \left( 2 \left( \alpha +\theta ^*+\nu \right) \left( 4 \alpha +4 \theta ^*+3 \root 3 \of {2} \root 3 \of {\lambda } \nu ^{2/3}\right) -8 \root 3 \of {2} \root 3 \of {\lambda } \nu ^{2/3} \left( \alpha +\theta ^*\right) \right) \left( \lambda \left( \frac{\nu ^2}{\left( \alpha +\theta ^*+\nu \right) ^2}-1\right) +\alpha +\theta ^*\right) ^2 \Bigg ] \Bigg \}. \end{aligned}$$

The expression is easy to evaluate numerically. Taking \(\lambda =1\) and \(\nu =4\) gives

$$\begin{aligned} \tilde{\xi }(\alpha ,\beta )=I_1(\alpha , \beta )+I_2(\alpha , \beta ), \end{aligned}$$

where

$$\begin{aligned} I_1(\alpha ,\beta )=\frac{1}{\left( -\beta ^3+3 \left( 2^{2/3}-4\right) \beta ^2+24 \left( 2^{2/3}-2\right) \beta +48\ 2^{2/3}-80\right) }\\ \times \frac{\left( 186+25 \root 3 \of {2}\right) \left( \alpha +2\ 2^{2/3}\right) \left( \alpha +2\ 2^{2/3}+8\right) (\beta +4)^2}{162\ 2^{2/3} \left( 3\ 2^{2/3}-5\right) ^7 \left( -\alpha ^2-4\ 2^{2/3} \alpha +\alpha +2\ 2^{2/3}-8 \root 3 \of {2}+4\right) ^2} \end{aligned}$$

and

$$\begin{aligned}&I_2(\alpha ,\beta )=\frac{\left( \left( \beta -3\ 2^{2/3}+4\right) (\beta +4)^2+16\right) ^{-2}}{324 \left( 5-3*2^{2/3}\right) ^6 \left( \alpha +2^{5/3}\right) \left( \left( \alpha +2^{5/3}-5\right) \left( \alpha +2^{5/3}\right) ^2+16\right) ^4 } \\&\quad \times \Bigg \{\bigg [5\ 2^{2/3} \alpha ^{11}+\left( 31 {2}^{7/3}-5^2\ 2^{5/3}\right) \alpha ^{10}+\left( 848-95*{2}^{13/3}-115\ 2^{2/3}\right) \alpha ^9\\&\qquad +4 \left( -4680+1269 \root 3 \of {2}+56*2^{2/3}\right) \alpha ^8-16 \left( -11412+4400 \root 3 \of {2}+4185*2^{2/3}\right) \alpha ^7\\&\qquad +64 \left( -24304-1215 \root 3 \of {2}+17476* 2^{2/3}\right) \alpha ^6-768 \left( -4875-8765 \root 3 \of {2}+9871*2^{2/3}\right) \alpha ^5\\&\qquad +2^9 \left( 40674-75618 \root 3 \of {2}+33715\ 2^{2/3}\right) \alpha ^4+2^{11} \left( -54132+32065 \root 3 \of {2}+8540\ 2^{2/3}\right) \alpha ^3\\&\qquad - 3*2^{13} \left( -4625-1418 \root 3 \of {2}+4002\ 2^{2/3}\right) \alpha ^2+3*2^{14} \left( 2155-2742 \root 3 \of {2}+805\ 2^{2/3}\right) \alpha \\&\qquad +2^{15} \left( -3780+450 \root 3 \of {2}+2029\ 2^{2/3}\right) \bigg ] (\beta +4)^2\left( -\beta ^3+3 \left( 2^{2/3}-4\right) \beta ^2+24 \left( 2^{2/3}-2\right) \beta +48\ 2^{2/3}-80\right) \\&\qquad -108 \root 3 \of {2} \left( \alpha +2^{5/3}-4\right) ^4 \left( \alpha +2^{5/3}\right) ^2 \left( \alpha +2^{5/3}+8\right) \left( -\alpha ^2-(2^{8/3}-1) \alpha +2^{5/3}- 2^{10/3}+4\right) ^2 (\beta +4)^4\Bigg \}. \end{aligned}$$

Observe that the Laplace transforms \(\tilde{\mu }(\alpha ,\beta )\) and \(\tilde{\xi }(\alpha ,\beta )\) can be inverted as they can be written as linear combinations of powers of \(1/(\mathrm{const}_\alpha +\alpha )\) and \(1/(\mathrm{const}_\beta +\beta )\) for various constants \(\mathrm{const}_\alpha \) and \(\mathrm{const}_\beta \).

Example 9

In this case, \(X(t)= B(t)-t,\) where B(t) is a standard Brownian motion. Note that this process is spectrally positive and spectrally negative. We apply the spectrally positive results. It is not hard to check that

$$\begin{aligned}\hat{\psi }({\vartheta })={\vartheta }+\frac{{\vartheta }^2}{2},\end{aligned}$$

so that \(\vartheta ^*= -1\) and \(\zeta ^*=-1/2.\) Furthermore, assumptions (SP) are satisfied and

$$\begin{aligned} \widehat{\Phi }(s)=-\left( 1+\sqrt{1+2s}\right) . \end{aligned}$$

It is a matter of straightforward computations now to obtain that

$$\begin{aligned}C_1(\alpha ,\beta )=-\frac{4 \sqrt{2}}{(\alpha +1)^2 (\beta +1)^2}\end{aligned}$$

and thus

$$\begin{aligned}\tilde{\mu }(\alpha ,\beta )=\frac{C_1(\alpha ,\beta )}{C_1(0,0)}=\left( \frac{1}{1 +\alpha }\right) ^2\left( \frac{1}{1+\beta }\right) ^2.\end{aligned}$$

In conclusion, the quasi-stationary distributions of Q(0) and Q(t) (conditioned that the busy period lasts longer than t, for large t) are both Erlang(2,1) with mean 2, whereas the stationary workload itself has an exponential distribution with mean 1/2; see [12, 21, 22]. Moreover, simple calculations lead to

$$\begin{aligned} C_3(\alpha ,\beta )=-\frac{8 \sqrt{2} (2 + \alpha (2 + \alpha ) + \beta (2 + \beta ))}{ (1 + \alpha )^4 (1 + \beta )^4}. \end{aligned}$$

Hence,

$$\begin{aligned}\tilde{\xi }(\alpha ,\beta )= \frac{3 \left( 2 \beta ^2+4 \beta +1\right) }{(\alpha +1)^2 (\beta +1)^4}-\frac{3}{(\alpha +1)^4 (\beta +1)^2}. \end{aligned}$$

One can observe that the second-order quasi-stationary measure \(\xi \) is given by

$$\begin{aligned}\xi (dx,dy)/dxdy=\frac{1}{2} x y e^{-(x+y)} \left( 12-x^2-y^2\right) . \end{aligned}$$

It is worth noting that \(\xi \) may be negative (as is the case here), since it pertains to the speed of convergence to the Yaglom limit.

6 Conclusions

In this paper we proved that the speed of convergence of Yaglom limits is very slow, namely of order 1/t. We also identified the prefactor in front of this speed, called the second-order quasi-stationary law. Additionally, we analysed in detail two examples.

Future research directions lie in Markov additive processes or affine processes. It is also interesting to analyse the Yaglom limit for multivariate Lévy processes where instead of the busy period one considers the first exit time from a positive quadrant in \(\mathbb {R}^d\) for \(d>1\).