\documentclass[12pt,reqno]{article}

\usepackage[usenames]{color}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{amscd}

\usepackage[colorlinks=true,
linkcolor=webgreen,
filecolor=webbrown,
citecolor=webgreen]{hyperref}

\definecolor{webgreen}{rgb}{0,.5,0}
\definecolor{webbrown}{rgb}{.6,0,0}

\usepackage{color}
\usepackage{fullpage}
\usepackage{float}

\usepackage{psfig}
\usepackage{graphics,amsmath,amssymb}
\usepackage{amsthm}
\usepackage{amsfonts}
\usepackage{latexsym}
\usepackage{epsf}

\usepackage{rotating}



%  Ward's macros begin
%
\newcommand{\bsq}{\vrule height .9ex width .8ex depth -.1ex}
\newcommand{\RR}{{\mathbb R}}
\newcommand{\ZZ}{{\mathbb Z}}
\newcommand{\sE}{{\cal E}}
\newcommand{\sA}{{\cal A}}
\newcommand{\sC}{{\cal C}}
\newcommand{\sD}{{\cal D}}
\newcommand{\sB}{{\cal B}}
\newcommand{\sG}{{\cal G}}
\newcommand{\sP}{{\cal P}}
\newcommand{\sQ}{{\cal Q}}
\newcommand{\sO}{{\cal O}}
\newcommand{\sU}{{\cal U}}
\newcommand{\sN}{{\cal N}}
\newcommand{\sW}{{\cal W}}
\newcommand{\sI}{{\cal I}}
\newcommand{\sR}{{\cal R}}
\newcommand{\sZ}{{\cal Z}}
\newcommand{\sS}{{\cal S}}
\newcommand{\sL}{{\cal L}}
\newcommand{\sM}{{\cal M}}
\newcommand{\sF}{{\cal F}}
\newcommand{\sT}{{\cal T}}
\newcommand{\sX}{{\cal X}}
\newcommand{\hsp}{\hspace*{\parindent}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\lf}{\lfloor}
\newcommand{\rf}{\rfloor}
\newcommand{\lc}{\lceil}
\newcommand{\rc}{\rceil}
\newcommand{\ra}{\rightarrow}
\newcommand{\deq}{\stackrel{\rm d}{=}}
\newcommand{\beql}[1]{\begin{equation}\label{#1}}
\newcommand{\eqn}[1]{(\ref{#1})}
\newcommand{\beq}{\begin{displaymath}}
\newcommand{\eeqno}{\end{displaymath}}
\newcommand{\vsp}{\vspace{.1in}}
\newcommand{\dd}{\ldots}
\newcommand{\af}{\alpha}
\newcommand{\la}{\lambda}
\newcommand{\La}{\Lambda}
\newcommand{\D}{\mathbb{D}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\bA}{{\mathbf A}}
\newcommand{\bzero}{{\mathbf 0}}
\newcommand{\bX}{{\mathbf X}}
\newcommand{\bY}{{\mathbf Y}}
\newcommand{\bD}{{\mathbf D}}
\newcommand{\bL}{{\mathbf L}}
\newcommand{\bS}{{\mathbf S}}
\newcommand{\bs}{{\mathbf s}}
\newcommand{\bN}{{\mathbf N}}
\newcommand{\bn}{{\mathbf n}}
\newcommand{\bC}{{\mathbf C}}
\newcommand{\bu}{{\mathbf u}}
\newcommand{\bc}{{\mathbf c}}
\newcommand{\bQ}{{\mathbf Q}}
\newcommand{\bB}{{\mathbf B}}
\newcommand{\bR}{{\mathbf R}}
\newcommand{\bV}{{\mathbf V}}
\newcommand{\bW}{{\mathbf W}}
\newcommand{\bI}{{\mathbf I}}
\newcommand{\bP}{{\mathbf P}}
\newcommand{\bfe}{{\bf e}}
\newcommand{\bK}{{\mathbf K}}
\newcommand{\bM}{{\mathbf M}}
\newcommand{\bU}{{\mathbf U}}
\newcommand{\bT}{{\mathbf T}}
\newcommand{\ua}{\uparrow}
\newcommand{\bZ}{{\mathbf Z}}
\newcommand{\bz}{{\mathbf z}}
\newcommand{\by}{{\mathbf y}}
\newcommand{\bw}{{\mathbf w}}
\newcommand{\bq}{{\mathbf q}}
\newcommand{\bb}{{\mathbf b}}
\newcommand{\ba}{{\mathbf a}}
\newcommand{\bm}{{\mathbf m}}
\newcommand{\bx}{{\mathbf x}}
\newcommand{\be}{{\mathbf e}}
\newcommand{\bone}{{\mathbf 1}}
\newcommand{\hv}{\hat{v}}
\newcommand{\hV}{\hat{V}}
\newcommand{\hf}{\hat{f}}
\newcommand{\hg}{\hat{g}}
\newcommand{\qandq}{\quad\mbox{and}\quad}
\newcommand{\qifq}{\quad\mbox{if}\quad}
\newcommand{\qforq}{\quad\mbox{for}\quad}
\newcommand{\qorq}{\quad\mbox{or}\quad}
\newcommand{\qasq}{\quad\mbox{as}\quad}
\newcommand{\qinq}{\quad\mbox{in}\quad}
\newcommand{\qforallq}{\quad\mbox{for all}\quad}
\newcommand{\df}{\displaystyle\frac}

%
%  Ward's macros end
%


\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{.1in}
\setlength{\evensidemargin}{.1in}
\setlength{\topmargin}{-.5in}
\setlength{\textheight}{8.9in}

\newcommand{\seqnum}[1]{\href{http://www.research.att.com/cgi-bin/access.cgi/as/~njas/sequences/eisA.cgi?Anum=#1}{\underline{#1}}}

\begin{document}

\begin{center}
\epsfxsize=4in
\leavevmode\epsffile{logo129.eps}
\end{center}

\begin{center}
\vskip 1cm{\LARGE\bf Integer Sequences from Queueing Theory}
\vskip 1cm \large
Joseph Abate\\
900 Hammond Road\\
Ridgewood, NJ 07450-2908\\
USA\\
%\href{mailto:mike@hempnallhouse.co.uk}{\tt mike@hempnallhouse.co.uk}\\
\ \\
Ward Whitt\\
Department of Industrial Engineering and Operations Research\\
Columbia University\\
New York, NY 10027-6699\\
USA\\
\href{mailto:ww2040@columbia.edu}{\tt ww2040@columbia.edu} \\
\end{center}

\vskip .2 in

\begin{abstract}
%Moment sequences associated with probability distributions can be
%sources of integer sequences.
  Operators on probability distributions can be expressed
  as operators on the associated moment sequences, and so
correspond to operators on integer sequences. Thus, there is an
opportunity to apply each theory to the other.
  Moreover, probability models can be sources
of integer sequences, both classical and new, as we show by
considering the classical $M/G/1$ single-server queueing model. We
identify moment sequences that are integer sequences.  We
establish connections between the $M/M/1$ busy period distribution and the Catalan and Schroeder numbers.
\end{abstract}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{defin}[theorem]{Definition}
\newenvironment{definition}{\begin{defin}\normalfont\quad}{\end{defin}}
\newtheorem{examp}[theorem]{Example}
\newenvironment{example}{\begin{examp}\normalfont\quad}{\end{examp}}
\newtheorem{rema}[theorem]{Remark}
\newenvironment{remark}{\begin{rema}\normalfont\quad}{\end{rema}}

\providecommand{\e}[1]{\ensuremath{\times 10^{#1}}}

\section{Introduction}\label{secIntro}
Given the well established link between integer sequences and combinatorics, it is surprising that
there are so few connections between integer sequences and queueing theory,
because two prominent scholars, John Riordan (1902--1988) and Lajos Tak\'{a}cs (1924-- \, ), were active in both
combinatorics and queueing theory \cite{R58,R62,R68,T62,T67,T90}.
Of course, these authors saw connections; e.g.,
the index of Riordan's queueing book \cite{R62} contains entries for Bell polynomials, binomial moments,
Catalan numbers, cumulant generating function, Lagrange expansion, and rooted trees.

In the OEIS \cite{S10}, a search on ``Riordan'' yields 1221 entries, while a search on ``queueing'' produces 2.
We think Riordan would want less disparity.
One exception is the entry submitted by A. Harel (\seqnum{A122525}) making connection to the Erlang delay formula associated with the classical $M/M/s$ queue.
A search on ``queue'' produces 15 entries, but these primarily focus on combinatorial problems \cite{A04,S05}
or a queue as a tool in computer science rather than queueing theory \cite{AB94};
queueing theory concerns stochastic process models describing congestion, e.g., the probability distribution of customer waiting times \cite{R62,T62}.

The purpose of this paper is to point out connections between the theories of probability and integer sequences,
and to exhibit some integer sequences that arise in queueing theory.
As a branch of probability theory, queueing theory has exploited
the classical analytical approach to probability theory
using transforms \cite{R62,T62}, whose series expansions involve the integer moments of the probability distributions
and close relatives.  The entire probability distribution
is characterized by the sequence of moments in great generality \cite{A65,DS97}.
Thus we were motivated to introduce an operational calculus for manipulating probability distributions on the positive halfline by either manipulating
the associated Laplace transforms or
the associated moment sequences \cite{AW96}; a quick overview is provided by \cite[Tables 1-3]{AW96}.
A more recent paper in the same spirit is \cite{FG00}.

The operational calculus for probability distributions on the positive halfline in the framework of moment sequences is closely related to operators
commonly used to analyze integer sequences.
%In the OEIS \cite{S10}, a search on ``moments'' yields only three consecutive related entries to
%the moments of a probability distribution (\seqnum{A0932526}).
Since many operators on integer
sequences can be applied to moment sequences arising in
probability theory, both when they are integers and when they are not, there is an opportunity for experts on
integer sequences to contribute to probability theory through moment sequences.
(The connection is also discussed in \cite{B09}.)  The
probability connection also provides concrete models where integer
sequences arise.

Here is how the present paper is organized.  We start in \S
\ref{secMoment} by reviewing moment sequences in probability
theory and relating operators on probability distributions to
operators on sequences.
Then in \S \ref{secMM1} we review the
classical $M/G/1$ single-server queueing model and identify
integer sequences arising in that context. We consider the $M/M/1$
busy period distribution, its stationary excess distribution and the equilibrium time to emptiness. In
that context we identify random quantities whose probability
density functions have the Catalan and Large Schroeder numbers as
moments. Finally, we consider the moment sequence for the $M/G/1$
steady-state waiting time. We give the proofs of Theorems \ref{th1} and \ref{th2} in
\S \ref{secProofs}.

\section{Moment Sequences in Probability}\label{secMoment}

  Let $Z$ be a nonnegative random variable with {\em cumulative distribution function} (cdf) $F$
and {\em probability density function} (pdf) $f$, i.e.,
\beq
F(t) \equiv P(Z \le t) = \int_{0}^{t} f(u) \, du, \quad t \ge 0,
\eeqno
where $\equiv$ means ``defined as.''
Let $\hat{f} (s)$ be the {\em Laplace transform} (LT) of $f$ (and thus $Z$) and let $\phi (x)$ be the associated {\em moment generating function} (mgf)
of $f$, i.e.,
\begin{eqnarray}\label{a2}
\hat{f}(s) & \equiv & \int_{0}^{\infty} e^{-st} f(t) \, dt \equiv E[e^{-sZ}] \qandq\nonumber \\
\phi (x)   & \equiv & \hat{f} (-x) = \int_{0}^{\infty} e^{x t} f(t) \, dt \equiv E[e^{xZ}],
\end{eqnarray}
where $s$ in \eqn{a2} is understood to be a complex number with positive real part, while $x$ in \eqn{a2}
is a positive real number.  The LT is always well defined; we assume that there exists $x^* > 0$ such that
$\phi (x) < \infty$ for $x < x^*$.

We obtain sequences by considering series expansions for the mgf $\phi$; in particular, we can write
\beql{a3}
\phi(x) = \sum_{n = 0}^{\infty} \frac{m_n x^n}{n!} = \sum_{n = 0}^{\infty} \mu_n x^n,
\eeq
where of course we must have $\mu_n = m_n/n!$, $n \ge 0$.
From the probability perspective, %are most interested in $m_n$, which is the $n^{\rm th}$ moment of the cdf $F$.
the object of primary interest is the cdf $F$, but the moments $m_n$ provide a useful partial characterization.
The moment sequence can be calculated and employed to derive other quantities of interest in probability models
\cite{ACLW95,AW96,CL96}.  There is a developing operational calculus for manipulating probability distributions via their moment sequences \cite{AW96}.

From the sequence perspective, we may instead regard the sequences
$\{m_n: n \ge 1\}$ and $\{\mu_n: n \ge 1\}$ as the objects of
primary interest.  The two perspectives meet with the mgf $\phi$.
From the sequence perspective, $\phi (x)$ arises as the {\em
generating function} (gf) of the sequence $\{\mu_n: n \ge 1\}$,
and we speak of $\{\mu_n\}$ as being its {\em coefficients}, while
$\phi (x)$ arises as the {\em exponential generating function}
(egf) of the sequence $\{m_n: n \ge 1\}$.  Of course, $\{m_n: n
\ge 1\}$ is always an integer sequence whenever $\{\mu_n: n \ge
1\}$ is, but not necessarily conversely. We will be giving
examples where both are integer sequences.
%In the theory of
%integer sequences, attention usually focuses on the gf and
%its coefficients $\{\mu_n\}$, but the egf is also used.

A simple example from probability theory is the exponential distribution with mean $m$, specified by
\beq
F(t) \equiv 1 - e^{-t/m}, \quad f(t) \equiv (1/m) e^{-t/m}, \quad t \ge 0, \qandq \phi (x) \equiv (1 - m x)^{-1},
\eeqno
which has associated sequences
\beq
m_n \equiv n! m^n \qandq \mu_n \equiv m^n, \quad n \ge 0.
\eeqno
It is immediate that $\{\mu_n\}$ and $\{m_n\}$ are both elementary integer sequences whenever $m$ is an integer.
Since the mean (first moment) is probabilistically only a scale parameter, depending on how the measurement units are defined,
it is natural to follow the convention that the first moment is $1$; here that gives $m_1 = 1$.  Then we obtain the fundamental integer sequences
$m_n = n!$ and $\mu_n = 1$, $n \ge 1$.

As a consequence of the relations outlined above, we see that results about probability distributions can be translated into
results about integer sequences, provided that the sequence $\{\mu_n\}$ or $\{m_n\}$ is indeed an integer sequence.  Similarly,
results about integer sequences can be translated into
results about probability distributions, provided that the egf or gf of the integer sequence is indeed the mgf of a
bonafide pdf.

We first observe that there is a natural probabilistic setting in which, not only is $\{m_n\}$ a moment sequence, but so is the
associated sequence $\{\mu_n\}$.  That occurs in spectral representations, where one pdf serves as a mixing pdf for another.
In particular, suppose that the pdf $f$ can be represented as a continuous mixture of exponential pdf's via
\beql{a31}
f(t) = \int_{\tau_1}^{\tau_2} y^{-1} e^{-t/y} g(y) \, dy, \quad t \ge 0,
\eeq
in which case we call $g$ the {\em mixing pdf} for $f$ \cite{AW88b}.
Append a superscript $f$ to $\{m_n\}$ and $\{\mu_n\}$ to denote dependence upon $f$.
We observe that the sequence $\{\mu_n\}$ is itself the moment sequence of the associated mixing pdf $g$.  Hence we call $\{\mu_n\}$ the {\em mixing moments} of $f$.
(We omit the elementary proofs in this section.)
\begin{proposition}{$($mixing moments$)$}\label{propMixing} For pdf's $f$ and $g$ related via {\em \eqn{a31}},
$m^g_n = \mu^f_n$, $n \ge 1$.
\end{proposition}

A canonical operation in probability theory is convolution.  If
two independent nonnegative random variables $Z_1$ and $Z_2$ with
pdf's $f_{Z{_1}}$ and $f_{Z{_2}}$ are added, then the sum $Z_1 +
Z_2$ has a pdf that is the {\em convolution} of the two component
pdf's, i.e.,
\beq
 f_{Z{_1} + Z{_2}} (t) = \int_{0}^{t}
f_{Z{_1}} (y) f_{Z{_2}} (t - y) \, dy, \quad t \ge 0.
\eeqno
One reason
transforms are so frequently used in probability is that they
convert convolution into a simple product; i.e., the associated
mgf's are related by
\beq
 \phi_{Z{_1} + Z{_2}} (x) =
\phi_{Z{_1}} (x) \phi_{Z{_2}} (x).
\eeqno
Because of the assumed stochastic independence, the moments are related by
the binomial theorem,
\beq
 m^{Z_{_1} + Z{_2}}_n =
\sum_{k=1}^{n} \binom{n}{k} m^{Z{_1}}_k m^{Z_{_2}}_{n-k}.
\eeqno
Thus, if $\{ m^{Z{_1}}_n\}$ and $\{ m^{Z{_2}}_n\}$ are integer
sequences, then so is $\{m^{Z{_1} + Z{_2}}_n\}$. For example, by
above, that occurs when $Z_1$ and $Z_2$ have exponential
or deterministic distributions with integer
means.


With this background, we can interpret \cite[Tables 1-3]{AW96}, which show various operators mapping one probability distribution
into another.  There are four columns:  The first column contains the name and notation for the operator; the second column shows how the operator acts on
LT's; the third column shows how it acts on pdf's; and the fourth column shows how it acts on moment sequences.
From the perspective of integer sequences, the fourth column shows how it acts on the coefficients of the egf;
we could then add a fifth column showing how it acts
on the corresponding coefficients of the gf.
Those familiar with integer sequences might want to translate the LT into the mgf by replacing $s$ with $-x$, and then interpret
that mgf as the egf of the given $k^{\rm th}$ moment.

In the rest of this section we highlight a few striking connections between operators on probability distributions, as in \cite{AW96},
and operators on integer sequences.
First, a standard operator on integer sequences is the simple {\em shift to the left}, e.g., converting
$1, 1, 2, 6, 22, 90, \ldots$ to $1, 2, 6, 22, 90, \ldots$.
  We now show that,
probabilistically, the simple shift applied to the coefficients $\mu_n$ corresponds to constructing the stationary-excess cdf
of a cdf on the positive halfline having mean $1$.

Given a nonnegative real-valued random variable $Z$ with cdf $F$ having
finite moments $m_k$, $k \ge 1$, let $Z_e$ be a random variable with the associated {\em stationary-excess cdf} $F_e$
(a.k.a. the equilibrium excess or stationary residual-life cdf),
defined by
\beql{a21}
F_e (t) \equiv P(Z_e \le t) \equiv \frac{1}{m_1}\int_{0}^{t} (1 - F(u)) \, du, \quad t \ge 0.
\eeq
The stationary-excess cdf frequently arises in renewal theory; see  of
\cite[Examples 7.16,7.17, 7.23, 7.24]{R07} and \cite{W85}; it appears in \cite[Table 2]{AW96}.
(A search on ``renewal theory'' in the OEIS gives two unrelated entries.)
For us, the important fact is that the
random variable $Z_e$ has moments
\beql{a22}
m_{e,k} \equiv E[Z_e^k] = \frac{m_{k+1}}{(k+1) m_1}, \quad k \ge 1.
\eeq
Hence, the transformation from a cdf of a nonnegative random variable to its associated
stationary-excess cdf produces a simple shift on the gf coefficients.
\begin{proposition}{$($the stationary-excess operator$)$}\label{propShift} Let $F$ be the cdf of a
nonnegative random variable with mean $1$, mgf $\phi$ in {\em \eqn{a2}} and
associated sequence of mixing moments $\{\mu_n\}$ in {\em \eqn{a3}} $($coefficients of $\phi (x)$ when it is regarded as a gf$)$.
The associated stationary-excess cdf $F_e$ in {\em \eqn{a21}} has mgf
$\phi_{e} (x) = (\phi (x) - 1)/x$.
The mixing moments of $\phi_e (x)$ $($coefficients of $\phi_e$ when it is regarded as a gf$)$ are
$\mu_{e,k} = \mu_{k+1}$, $k \ge 1$.
\end{proposition}

A similar relationship holds for the {\em stationary-lifetime operator}, mapping a pdf $f$ into the associated pdf
\beql{a25}
f_s (t) \equiv \frac{t f(t)}{m_1}, \quad t \ge 0.
\eeq
The stationary-lifetime pdf also frequently arises in renewal theory, e.g., \cite[\S 7.7]{R07}, and also appears in \cite[Table 2]{AW96}.
For us, the important fact is that the
moments are related by
\beq
m_{s,k} \equiv \frac{m_{k+1}}{m_1}, \quad k \ge 1,
\eeqno
just like \eqn{a22} without the $k+1$ in the denominator.
Hence, the transformation from a pdf of a nonnegative random variable to its associated
stationary-lifetime pdf produces a simple shift on the moments $m_n$ (the egf coefficients).
\begin{proposition}{$($the stationary-lifetime operator$)$}\label{propShift2}
Let $f$ be the pdf of a nonnegative random variable with mean $1$, mgf $\phi$ in {\em \eqn{a2}} and
associated sequence of moments $\{m_n\}$ in {\em \eqn{a3}}.  The associated stationary-lifetime pdf $f_s$ in {\em \eqn{a25}} has mgf
$\phi_{s} (x) = \phi' (x)$.
The coefficients of $\phi_s$ when it is regarded as an egf are
$m_{s,k} = m_{k+1}$, $k \ge 1$.
\end{proposition}

We now turn to another basic probability operator, which is conveniently related to a
continued fraction representation of the mgf \cite{AW99a,AW99b, B09, CBVWJ08}. % and thus to Hankel transforms; see \S \ref{secHankel} below.
If $\hat{f}$ is the LT of a pdf $f$, then the associated {\em exponential mixture pdf} has LT
\beql{a41}
\hat{f}_{\sE\sM} (s) \equiv (1 + s \hat{f}(s))^{-1};
\eeq
it appears in \cite[Table 3]{AW96}.
The special case of exponential mixtures of inverse Gaussian (EMIG) distributions is discussed in \cite[\S 8]{AW96} and in \cite{AW99b}.

The probabilistic exponential mixing operator has a simple manifestation in the continued fraction representation of the mgf,
%and the Hankel transform of the coefficients,
regarding that mgf as a gf.
Starting with the LT $\hat{f} (s)$ of a pdf $f$, if we represent the associated mgf as a formal power series
by
\beql{power}
\phi (x) \equiv \hat{f} (-x) = 1 + \mu_1 x + \mu_2 x^2 + \mu_3 x^3 + \mu_4 x^4 + \mu_5 x^5 + \ldots ,
\eeq
then the {\em corresponding continued fraction} (CF) is
\beql{CF}
\phi (x) = \frac{1}{1-} \frac{h_1 x}{1-} \frac{h_2 x}{1-} \frac{h_3 x}{1-} \frac{h_4 x}{1-} \ldots
\eeq
where $h_1 \equiv \mu_1$, $h_2 \equiv (\mu_2 - \mu_1^2)/\mu_1$, etc.  When $h_n > 0$ for all $n$
we have an $S$-fraction and the underlying pdf $f$ is {\em completely monotone} (CM) \cite{AW99a}.

\begin{proposition}{$($the exponential-mixture operator$)$}\label{propShift3} Let $f$ be the pdf of a nonnegative random variable with mean $1$,
mgf $\phi$ in {\em \eqn{a2}} with associated sequence of CF coefficients $\{h_n\}$ in {\em \eqn{CF}}.
The associated exponential-mixture pdf $f_{\sE\sM}$ with LT in {\em \eqn{a41}} has CF coefficients
\beq
h_{\sE\sM,1} = 1, \quad h_{\sE\sM,n} = h_{n-1}, \quad n \ge 2;
\eeqno
i.e., it produces a shift to the right.  The inverse exponential mixture operator
{\em \cite[(7.3), p. \ 94]{AW96}} gives the corresponding shift to the left.
\end{proposition}

\paragraph{Proof.}  From the transform expression in \eqn{a41}, the conclusion is immediate:  Given the CF
representation for $\hat{f}(s)$, it is immediate that the corresponding CF for $(1 + s \hat{f}(s))^{-1}$
shifts the coefficients one to the right; i.e., we write the CF for $(1 + s \hat{f}(s))^{-1}$ as
$1/(1+s \hat{f} (s))$, inserting the CF for $\hat{f} (s)$.~~~\bsq


%In the theory of integer sequences, the operator corresponding to the shift in the coefficients of the CF is well known.
%In Proposition \ref{propShift3} we are giving a probabilistic counterpart.




\section{Queueing Examples}\label{secMM1}

\subsection{The $M/G/1$ Model}

The $M/G/1$ queue is a basic model in queueing theory, usually discussed in queueing textbooks; e.g., \cite[\S 8.5]{R07}.  There is a single server
with unlimited waiting room.   Customers arrive according to a Poisson process (the first $M$, for Markov) with rate $\lambda$, $0 < \lambda < \infty$.
If the system is empty, then the customer goes immediately into service; otherwise the customer waits in queue.
The successive service times come from a sequence of independent and identically distributed (i.i.d.) random variables with cdf $G$
having mean $1/\mu$, $0 < \mu < \infty$.  We will mostly consider the easiest case, in which the service-time cdf $G$ is exponential;
then the model is denoted by $M/M/1$.

A waiting customer enters service immediately upon service completion.
Let $Q(t)$ be the number of customers in the system at time $t$ for $t \ge 0$.  In the $M/M/1$ model, the stochastic process $Q \equiv \{Q(t): t \ge 0\}$ is
a birth-and-death stochastic process, with constant birth rate $\lambda$ and constant death rate $\mu$.  Let $\rho \equiv \lambda/\mu$ be the {\em traffic intensity}.
If $\rho < 1$, then $P(Q(t) = j|Q(0) = i) \ra (1 - \rho)\rho^j$ as $t \ra \infty$ for each $i$ and $j$; i.e.,
$Q(t)$ converges in distribution to a geometric distribution on the nonnegative integers, having mean $\rho/(1- \rho)$.
We assume that $\rho < 1$, under which the system is said to be stable.
(If $\rho \ge 1$, then $P(Q(t) \le j|Q(0) = i) \ra 0$ as $t \ra \infty$ for each $i$ and $j$.)
For the $M/G/1$ model, the steady-state distribution of $Q(t)$ is characterized by the Pollaczek-Khintchine transform \cite{R62,T62}.

\subsection{The $M/M/1$ Busy Period}

A busy period is the time from the arrival of a customer finding an empty system until the system is empty again \cite[\S 4.8]{R62}.  The first passage time of $Q$
from any state $j> 0$ to $j - 1$ is distributed as a busy period.
Without loss of generality, we can measure time in units of mean service times, so that we let $\mu \equiv 1$.  Then the model has the single parameter $\rho \, (=\lambda)$.
Let $X_{\rho}$ denote a random variable with the busy period distribution, as a function of the traffic intensity $\rho$.  Let $B_{\rho}$ be the
cdf of $X_{\rho}$, i.e., $B_{\rho} (t) \equiv P(X_{\rho} \le t)$, $t \ge 0$,
and let $b_{\rho}$ be the associated pdf.  It turns out that
\beq
b_{\rho} (t) \equiv \frac{1}{t\sqrt{\rho}}e^{-(1 + \rho)t} I_{1}(2t\sqrt{\rho}), \quad t \ge 0,
\eeqno
where $I_{1} (t)$ is the Bessel function of the first kind \cite[(39), p. 63]{R62}.
The LT of $b_{\rho}$ (and thus $X_{\rho}$) is
\beq
\hat{b}_{\rho} (s) \equiv \int_{0}^{\infty} e^{-st} b_{\rho} (t) \, dt \equiv E[e^{-sX_{\rho}}] = \frac{1 + \rho + s - \sqrt{(1 + \rho +s)^2 - 4 \rho}}{2\rho},
\eeqno
\cite[(38), p. 63]{R62}.  As usual, the moments can be obtained by differentiating the Laplace transform.  The first two moments are
$E[X_{\rho}] = (1 - \rho)^{-1}$ and $E[X_{\rho}^2] = 2(1 - \rho)^{-3}$.

\subsection{The Busy-Period Moment Sequence After Scaling}

Our goal here is to obtain interesting integer sequences from the sequence of successive integer moments of a busy period $X_{\rho}$.
To obtain integer sequences, we first perform a change of variables,
introducing $\sigma \equiv \rho/(1 - \rho)$ (the mean steady state number in system) or, equivalently, $\rho \equiv \sigma/(1 + \sigma)$.
Notice that the first two moments become $1 + \sigma$ and $2 (1 + \sigma)^3$.  It turns out that the entire moment sequence becomes an integer sequence
whenever $\sigma$ is an integer.  Hence, from this first step, we obtain an entire sequence of integer sequences.

To seek simple integer sequences, we further scale these moment sequences all to have mean (first moment) $1$.
We remark that this final spatial scaling also plays a role in understanding the performance of
the $M/M/1$ queue as the traffic intensity $\rho$ increases toward its critical value $1$.  With appropriate scaling
of both time and space, the stochastic process $Q$ approaches reflected Brownian motion
with negative drift, while the busy period distribution has interesting behavior, in which both small values and large values play a role
\cite{AW88a,AW95,W02}.

In terms of random variables, let
\beql{b3}
Y_{\sigma} \equiv \frac{X_{\sigma/(1+\sigma)}}{1 + \sigma} =
 (1 - \rho) X_{\rho}, \quad \mbox{where} \quad \rho \equiv \frac{\sigma}{1 + \sigma} \qandq \sigma = \frac{\rho}{1 - \rho}.
\eeq
Let $b(x; \sigma)$ be the mgf of $Y_{\sigma}$.  From above,
\beql{b4}
b(x; \sigma) \equiv E[e^{x Y_{\sigma}}] = \frac{1 + 2 \sigma - x - \Psi (x)}{2 \sigma}, \quad \mbox{where} \quad  \Psi (x) \equiv \sqrt{1 - 2(1 + 2 \sigma)x + x^2}.
\eeq
Let the associated moments of
$Y_{\sigma}$ be
$m_n (\sigma) \equiv E[Y_{\sigma}^n]$, $n \ge 1$,
where $m_1 (\sigma) = 1$ for all $\sigma > 0$.  These moments, divided by $n!$, are the coefficients of the series expansion of $b(x; \sigma)$
\beql{b6}
b(x; \sigma) = \sum_{n = 0}^{\infty} \frac{m_n (\sigma) x^n}{n!}
=  \sum_{n = 0}^{\infty} b_n (\sigma) x^n, \quad \mbox{where} \quad b_n (\sigma) \equiv \frac{m_{n} (\sigma)}{n!}.
\eeq
The coefficients $b_n (\sigma)$ (and thus also
the moments $m_n (\sigma)$) are polynomials in $\sigma$ (and thus integers when $\sigma$ is an integer); the first few are:
$b_0 (\sigma) = 1$, $b_1 (\sigma) = 1$, $b_2 (\sigma) = 1 + \sigma$ and $b_3 (\sigma) = 1 + 3 \sigma + 2 \sigma^2$.

More specifically, we now relate the coefficients $b_n (\sigma)$ in \eqn{b6} to the Catalan numbers, denoted by $C_k$, (\seqnum{A000108})
starting with $1, 1, 2, 5, 14, 42, 132, 429, 1430$; in particular,
\beql{Catalan}
C_n \equiv \frac{1}{n+1} \binom{2n}{n} \qandq C_{n+1} = \sum_{i=0}^{n} C_i C_{n-i}, \quad n \ge 1.
\eeq
  The Catalan numbers can be
characterized in terms of their generating function
\beql{b7}
c(y) \equiv \sum_{k=0}^{\infty} C_k y^k \equiv \frac{1 - \sqrt{1 - 4y}}{2y}.
\eeq
%The Hankel transforms of this sequence give $1, 1, 1, \ldots$.
\begin{theorem}{$($mixing moments of the busy period pdf$)$}\label{th1}  For $n \ge 1$,
\beql{b8}
b_{n+1} (\sigma) = \sum_{k = 0}^{n} \binom{n+k}{n-k} C_k \sigma^k,
\eeq
where $b_n (\sigma)$ is defined in {\em \eqn{b4}} and {\em \eqn{b6}}, and $C_k$ is the $k^{\rm th}$ Catalan number in {\em \eqn{Catalan}}.
In addition, there is a convenient recurrence relation, with
$b_{0} (\sigma) \equiv b_{1} (\sigma) \equiv 1$ and
\beql{b8a}
(n + 1) b_{n+1} (\sigma) = (2n -1)(1 + \sigma) b_n (\sigma) - (n-2) b_{n-1} (\sigma).
\eeq
\end{theorem}
We provide a proof in \S \ref{secProofs}.  For $\sigma = 1,2, 3$, the coefficient sequences are:
\begin{eqnarray}
\{b_n (1)\} & = & 1, \quad 1, \quad  2, \quad  6, \quad  22, \quad  90, \quad  394, \quad  (\seqnum{A155069}) \nonumber \\
\{b_n (2)\} & = & 1, \quad 1, \quad  3, \quad  15, \quad  93, \quad 645, \quad 4791, \quad  (\seqnum{A103210}) \nonumber \\
\{b_n (3)\} & = & 1, \quad 1, \quad  4, \quad  28, \quad  244, \quad  2380, \quad  24868, \quad  (\seqnum{A103211}) \nonumber
\end{eqnarray}
Sequence $\{b_n (1)\}$ (\seqnum{A155069}) is a relatively recent addition to the OEIS \cite{S10}, having the title ``Expansion of
$(3 - x - \sqrt{1 - 6x +x^2})/2$ in powers of $x$.''  We provide a model context.  However, sequence
$\{b_n (1)\}$ shifted one to the left, i.e., $1, 2, 6, 22, 90, \ldots$ is 
(\seqnum{A006318}), corresponding to the famous
``Large Schroeder numbers.''  The ``little Schroeder numbers'' in (\seqnum{A001003}) are obtained by dividing \seqnum{A006318} by $2$, i.e.,
the sequence $1, 1, 3, 11, 45, \ldots$.

The moment sequences themselves, $m_n (\sigma) \equiv n! b_n (\sigma)$ are of course
also integer sequences, but they are evidently not in the OEIS \cite{S10}.  For example, the sequences
with terms $n! b_n (1)$ and $n! b_n(1/2)$, $n \ge 0$, yielding $1,1,4,36, 528, 10800$  and $1,1,3, 18,171, 2250$, respectively, are not found.

The busy period is the first passage time of $Q$ from $1$ to $0$.  Since the first passage time pdf from each state $k$ to $0$ is the $k$-fold convolution
of the busy period pdf \cite[Theorem 3.1]{AW88a}), the moments and mixing moments of these pdf's also generate integer sequences for integer $\sigma$.

Since we find regularity by multiplying $\{b_n (\sigma)\}$ by $1/2 = \sigma/(1 + \sigma)$ when $\sigma = 1$, we are motivated to consider
the associated sequences $\{\sigma b_{n+1} (\sigma)/(1 + \sigma): n \ge 1\}$ for integers $\sigma > 1$.
\begin{corollary}\label{cor1} For $n \ge 1$,
\beq
\frac{\sigma b_{n+1} (\sigma)}{1 + \sigma} = \sum_{k = 0}^{n} (-1)^{n+k} \binom{n+k}{n-k} C_k (1 + \sigma)^k.
\eeqno
\end{corollary}

We have already observed that Corollary \ref{cor1} gives the little Schroeder numbers for $\sigma = 1$; it is also discussed in
\cite[Problem 15, p. 168]{R68}.

\subsection{The Busy-Period Stationary-Excess Distribution}\label{secExcess}

Let $B_e$ denote the busy period stationary-excess cdf associated with the busy-period cdf $B_{\rho}$, defined as in \eqn{a21}
after applying the scaling in \eqn{b3}.  Let $b_e$ be the associated stationary-excess pdf.
 We can apply Proposition \ref{propShift} to characterize the mixing moments $b_{e,n} (\sigma)$ of $b_e$.
\begin{corollary}{$($mixing moments of the busy-period stationary-excess pdf$)$}\label{cor2} The $M/M/1$ busy-period stationary excess pdf $b_e$
has mgf
\begin{eqnarray}\label{be}
b_e (x, \sigma) & \equiv &\int_{0}^{\infty} e^{x t} b_e (t) \, dt \equiv \sum_{n=0}^{\infty} b_{e,n} (\sigma)x^n = \frac{b(x; \sigma) - 1}{x} \nonumber \\
& = & \frac{2}{1 - x + \sqrt{1 - 2(1+ 2\sigma)x + x^2}} \nonumber \\
& = & \frac{c(\sigma x/(1-x)^2)}{1-x} = \frac{2}{1 - x + \Psi(x)},
\end{eqnarray}
where $c(y)$ is the generating function of the Catalan numbers in {\em \eqn{b7}} and
$\Psi$ is defined in {\em \eqn{b4}}.
For $n \ge 1$,
\beql{ex1}
b_{e,n} (\sigma) = b_{n+1} (\sigma), n\ge 1,
\eeq
 for which an expression is given in Theorem {\em \ref{th1}}.  The mean of $b_e$ is $\sigma$.
\end{corollary}


\subsection{Large Schroeder and Catalan Numbers as Moments of PDF's}

In this section we identify the pdf's whose moments are (i) the Large Schroeder numbers and (ii) the Catalan numbers.
By Corollary \ref{cor2} and our observation after Theorem \ref{th1}, the sequence $\{b_{e,n} (1)\}$ coincides with the Large Schroeder numbers (\seqnum{A006318}).
From \cite[Theorem 5.1 and Corollary 5.2.1]{AW88a} and \cite[Theorem 4.1]{AW88b}, we can obtain the spectral representation for the pdf $b_e$.
(We need to account for the different scaling of time used there.)  That identifies the desired pdf and, by Proposition \ref{propMixing}, the mixing moments.

\begin{corollary}{$($Large Schroeder numbers as moments$)$}\label{cor3} The large Schroeder numbers arise as the moments of the mixing pdf
associated with the stationary-excess busy-period pdf $b_e$ when $\sigma = 1$.
The mixing pdf for $b_e$ as a function of $\sigma$ is
 \beq
f(y; \sigma) = \frac{\sqrt{(\tau - y)(y - \tau^{-1})}}{2 \sigma \pi y}, \quad \frac{1}{\tau} < y < \tau,
\eeqno
where $\tau \equiv 1 + 2\sigma + 2\sqrt{\sigma(1 + \sigma)}$.
\end{corollary}

We introduced the Catalan numbers before Theorem \ref{th1}.  The Catalan numbers are related
to reflected Brownian motion with
drift $-1$ and diffusion coefficient $1$, denoted by $\{R(t): t\ge0\}$.
It is the limit of the $M/M/1$ queue length process as $\rho \uparrow 1$ with appropriate scaling of
time and space \cite{AW87,W02}.  Since $E[R(t)|R(0) = 0]$ is nondecreasing,
$H_1 (t) \equiv E[R(t)|R(0) = 0]/E[R(\infty)]$, $t \ge 0$, is a cdf.
By \cite[Corollary 1.5.2]{AW87}, \cite[\S 7]{AW96} and \cite[(8.7)]{AW99a},
the pdf $h_1$ of $H_1$ has LT $\hat{h}_1 (s)$, which can be characterized as the unique fixed point of the exponential mixture operator, i.e.,
\beql{h2}
\hat{h}_1 (s) = \frac{1}{1 + s \hat{h}_1 (s)}.
\eeq
Together with the spectral representation given in \cite[Theorem 4.1]{AW88b},
that implies the following result.  Without the probability model context, the result already appears in the commentary on
(\seqnum{A000108}) \cite{P02}.
\begin{theorem}{$($Catalan numbers as moments$)$}\label{thCat} The generating function $c(x)$ of the
Catalan numbers in {\em \eqn{Catalan}} and {\em \eqn{b7}} coincides
with the mgf of $h_1$, $\hat{h}_1 (-x)$, associated with the LT in {\em \eqn{h2}}.
As a consequence, the Catalan numbers arise as the moments of the mixing density of $h_1$,
\beql{h3}
f(y) = \frac{\sqrt{4 - y}}{2\pi \sqrt{y}}, \quad 0 < y < 4.
\eeq
\end{theorem}

If we take the Laplace transform of $f$ in \eqn{h3} and replace $s$ by $-x$, then we obtain the mgf,
which is the egf $\tilde{c} (x)$ of the Catalan numbers.
\begin{corollary}\label{corEGF}
The egf of the Catalan numbers is
\beq
\tilde{c} (x) \equiv \sum_{n=0}^{\infty} \frac{C_n x^n}{n!} = e^{2x}(I_{0} (2x) - I_{1} (2x)),
\eeqno
where $I_0$ and $I_1$ are Bessel functions.
\end{corollary}

\subsection{The $M/M/1$ Equilibrium Time to Emptiness}

From a sequence perspective, the shifting by $1$ and multiplying by $\sigma/(1+\sigma)$ following Theorem \ref{th1} and Corollary \ref{cor1}
lead us to consider the mgf
\beql{b10}
p(x; \sigma) \equiv \frac{1}{1 + \sigma} + \frac{\sigma}{1 + \sigma}\left(\frac{b(x; \sigma) - 1}{x}\right).
\eeq

The function $p(x; \sigma)$ in \eqn{b10} turns out to be the moment generating function of the equilibrium time to emptiness
in the $M/M/1$ queue, i.e., the first passage time to state $0$ by the stochastic process $Q$, assuming that $Q$ starts at time $0$ according to
its (geometric) steady-state distribution.
Hence, with probability $(1-\rho) = 1/(1 + \sigma)$ the process starts at $0$, so the first passage time is $0$.
As a consequence, there is an atom at $0$ with mass $1/(1 + \sigma)$.  Conditional on starting at a positive value, the equilibrium time to emptiness
coincides with the stationary-excess pdf $b_e$ of the busy-period pdf $b$ \cite[Theorem 3]{AW94}.  Additional characterizations appear
in \cite[Theorem 3.3]{AW88a}.  We make a connection to the exponential mixture operator in \eqn{a41}.
Paralleling \eqn{b6}, we let $p_n (\sigma)$ be the coefficient of $p(x; \sigma)$ as a gf, i.e., writing
$p (x;\sigma) \equiv \sum_{n=0}^{\infty} p_n (\sigma) x^n$.
%The following draws on the reversion of series, \cite{R68}, p. 149-151.
We provide a proof in \S \ref{secProofs}.

\begin{theorem}{$($equilibrium time to emptiness$)$}\label{th2}
The mgf $p(x;\sigma)$ in {\em \eqn{b10}} can be expressed as
\beql{b15}
p(x; \sigma)  =  \frac{c((1+\sigma)x/(1+x)^2)}{1+x} = \frac{2}{1 + x + \Psi(x)},
\eeq
where $c(y)$ is the generating function of the Catalan numbers in {\em \eqn{b7}} and
$\Psi$ is defined in {\em \eqn{b4}}.  Hence,
\beql{b15a}
p_n (\sigma) = \frac{\sigma b_{n+1}(\sigma)}{1+\sigma} =  \frac{\sigma}{1 + \sigma} \sum_{k=0}^{n} \binom{n+k}{n-k} C_k \sigma^k, \quad n \ge 0,
\eeq
and the following recursion can be applied with $p_0 (\sigma) \equiv 1$ and $p_1(\sigma) \equiv \sigma$
\beql{b15b}
(n+2) p_{n+1} (\sigma) = (2n+1)(1+2\sigma) p_n (\sigma) - (n-1) p_{n-1}(\sigma), \quad n \ge 1.
\eeq
In addition, $b(x; \sigma)$ can be expressed as an exponential mixture of $p(x; \sigma)$, i.e.,
\beql{b16}
b(x; \sigma) = \frac{1}{1 - x p(x; \sigma)}.
\eeq
\end{theorem}


From \eqn{b15a} and \eqn{b15b},
 the sequence $\{p_n (\sigma): n \ge 1\}$ is an integer sequence for each positive integer $\sigma$.
For $\sigma = 1,2, 3$, the coefficient sequences are:
\begin{eqnarray}
\{p_n (1)\} & = & 1, \quad 1, \quad  3, \quad  11, \quad  45, \quad  197, \quad  (\seqnum{A001003}) \nonumber \\
\{p_n (2)\} & = & 1, \quad 2, \quad  10, \quad  62, \quad  430, \quad 3194, \quad  (\seqnum{A107841}) \nonumber \\
\{p_n (3)\} & = & 1, \quad 3, \quad  21, \quad  183, \quad  1785, \quad 18651, \quad   (\seqnum{A131763}). \nonumber
\end{eqnarray}

We conclude this section by remarking that in the theory of integer sequences there also is a convenient {\em invert operator},
which can be expressed for LT's via
\begin{eqnarray}\label{invert}
\hat{g} (s) & \equiv & \sI (\hat{f}(s)) \equiv \mbox{Excess} (\mbox{exp mixture}(\hat{f} (s)) = \hat{f}_{\sE\sM, e} (s) \nonumber \\
            & = &\mbox{Excess} \left(\frac{1}{1 + s \hat{f}(s)}\right) = \frac{1}{s} \left(1 - \frac{1}{1 + s \hat{f} (s)}\right) \nonumber \\
            & = & \frac{\hat{f}(s)}{1 + s \hat{f} (s)}.
\end{eqnarray}

The following result links the three mgf's $b(x; \sigma)$, $b_e (x; \sigma)$ and $(p (x; \sigma)$; see Corollary 5.1 of \cite{AW88a}.
We also consider the excess of the excess, denoted by $b_{e.e}$ \cite{W85}.
\begin{corollary}\label{cor4}  The mgf's $b_e (x; \sigma)$, $p (x; \sigma)$ and $b(x; \sigma)$ are related by
\beq
b_e (x; \sigma) = \sI (p (x; \sigma)) = p(x; \sigma) b(x; \sigma)
\eeqno
for the invert operator in {\em \eqn{invert}} and
\beq
b_{e,e} (x; \sigma) = p_e (x; \sigma).
\eeqno
\end{corollary}


Directly from the final expression in \eqn{be}, we get a heavy-traffic limit for the scaled mgf
with our scaling in \eqn{b3}, with the limit
being the gf of the Catalan numbers.  Let $X_e (\sigma)$ be a random variable with mgf $b_e (x; \sigma)$;
i.e., $b_e (x;\sigma) \equiv E[e^{x X_e (\sigma)}]$.  Let $\Rightarrow$ denote convergence in distribution \cite{W02}.
\begin{corollary}\label{cor5}  As $\sigma \ra \infty$,
\beq
 b_e (x/\sigma; \sigma) \ra  \frac{2}{1 - \sqrt{1 - 4x}} = \frac{1 - \sqrt{1 - 4x}}{2x} \equiv c(x)
\eeqno
for $c$ in \eqn{b7}, so that
\beq
\frac{X_e (\sigma)}{\sigma} \Rightarrow X \qasq \sigma \ra \infty,
\eeqno
where $c(x) = E[e^{x X}]$.
\end{corollary}
A related heavy-traffic limit for the time-scaled busy-period stationary-excess pdf was obtained in
\cite[Theorem 1]{AW95}.  (They are consistent.)

\subsection{Continued Fractions and Hankel Transforms}

We now look at the three probability mgf's $b(x; \sigma)$, $b_e (x,; \sigma) \equiv (b(x;\sigma) - 1)/x$ and $p(x; \sigma)$ in
\eqn{b4}, \eqn{be} and \eqn{b10} from the perspective of continued fractions and Hankel transforms.
From the previous subsections, we know that the associated three random variables are closely related.
That is shown again through this new perspective.  First, the three mgf's can be represented as $S$ (Stieltjes) fractions.
(In \cite{AW99a} we found that $S$ fractions frequently arise in CF's associated with birth-and-death stochastic processes.)
The CF coefficients are given in Table \ref{tb1}.
\begin{table}[h!]
\begin{center}
\begin{tabular}{|r||llllll|}
\hline \hline
mgf             & $h_1$ & $h_2$ & $h_3$ & $h_4$ & $h_5$ & $h_6$ \\ \hline
$b(x; \sigma)$     & 1         & $\sigma$    & $1 + \sigma$ & $\sigma$     & $1+\sigma$ & $\sigma$ \\
$b_e(x; \sigma)$  & $1+ \sigma$ & $\sigma$    & $1 + \sigma$ & $\sigma$     & $1+\sigma$ & $\sigma$ \\
$p(x;\sigma)$     & $\sigma$    & $1+ \sigma$ & $\sigma$     & $1 + \sigma$ & $\sigma$   & $1+\sigma$ \\
\hline \hline
\end{tabular}
\caption{The coefficients of the continued fractions representing the three $M/M/1$ mgf's.
The pattern repeats in each row.}\label{tb1}
\end{center}
\end{table}

For the generating functions of our $M/M/1$ queueing examples, the determination of the CF
coefficients is simple because we can apply the algebraic identity
\beq
 \frac{\alpha }{1+} \frac{\gamma}{1+} \frac{\alpha}{1+} \frac{\gamma}{1+} \ldots
 = \frac{1}{2}\left(\sqrt{1 + 2(\gamma + \alpha) + (\gamma - \alpha)^2} - 1 - (\gamma - \alpha)\right)
\eeqno
for constants $\alpha$ and $\gamma$;
i.e., we can solve the equation
\beq
\zeta = \frac{\alpha}{1+} \frac{\gamma}{1+ \zeta}
\eeqno
for $\zeta$.

The Hankel transform of an integer sequence provides a useful partial characterization;
it is a many-to-one function mapping an integer sequence into another integer sequence.
(For an example of non-uniqueness, see Corollary \ref{cor6} below.)
For example, the Hankel transform of the Catalan numbers is the sequence $1,1,1, \dots$ \cite{B09, L01}.
Following \cite{CBVWJ08, DS97}, starting from
a sequence $\{\omega_n: n \ge 0 \} \equiv \omega_0, \omega_1, \omega_2, \omega_2, \dots$
with $\omega_0 \equiv 1$, let the
{\em Hankel matrix} $M^{(n)}$ be the $(n+1) \times (n+1)$ symmetric matrix with elements $M^{(n)}_{i,j} \equiv \omega_{i+j-2}$,
$0 \le i \le n$, $0 \le j \le n$.
(The first row contains the first $n+1$ elements and $M_{n+1,n+1} \equiv \omega_{2n}$.
Let $H_{2n} \equiv det(M^{(n)})$, the {\em even  Hankel determinant}.
Let the {\em Hankel transform} of the sequence $\{\omega_n: n \ge 0 \}$ above be the sequence
$\{H_{2n}: n \ge 0\}$; it starts with $H_0 = 1$.

One way to compute the Hankel transform of a sequence $\{\omega_n: n \ge 0 \}$ is to first
determine the corresponding continued fraction from the formal power series; i.e.,
starting from \eqn{power}, we determine \eqn{CF} above.  The coefficients $h_n$ appearing in \eqn{CF}
can be determined by applying the normalized Viskovatov algorithm, as given on \cite[p. 112]{CBVWJ08}.
Then we apply the iteration
\beq
H_{2n} = \left(\prod_{i = 1}^{2n} h_i \right) H_{2n - 2}, \quad n \ge 1;
\eeqno
see Theorem 1.4.10 on  of \cite[p. 23]{DS97}.  Another way to arrive at this result is via the even contraction of the CF in \eqn{CF},
as given in \cite[(12.3), p. 270]{B09}; note that $\beta_n$ in \cite{B09} is equal to $h_{2n-1} h_{2n}$, $n \ge 1$, in our notation.
Then (12.2) of \cite[(12.2)]{B09} and our iteration yield the same result.


\begin{corollary}\label{cor6}
The Hankel transforms of the sequences $\{b_n (\sigma)\}$ in Theorem {\em \ref{th1}} , $\{b_{e,n} (\sigma)\}$ in Corollary {\em \ref{cor2}}
and $\{p_n (\sigma\}$ in Theorem {\em \ref{th2}} are
\begin{eqnarray}
H_{2n} (b)   & = & \sigma^{\nu(n)}(1 + \sigma)^{\nu(n) - n}, \nonumber \\
H_{2n} (p)   & = & H_{2n} (b_e) = (\sigma + \sigma^2)^{\nu(n)}, \nonumber
\end{eqnarray}
where $\nu (n) \equiv n(n+1)/2$.
\end{corollary}

\subsection{The Stationary Waiting Time in the $M/G/1$ Queue}

The mgf $w(x; \sigma) \equiv E[e^{xW}]$ of the steady-state waiting time $W$ in the $M/G/1$ queue
can be expressed in terms of the mgf $g_e (x; \sigma)$ of the stationary-excess
of the general service-time distribution using the Pollazcek-Khintchine transform (again assuming that the mean service time is $1$) by
applying the random sum operator from \cite[Table 1]{AW96}, yielding
\beql{w1}
w(x; \sigma) = \frac{1-\rho}{1 -\rho g_e (x; \sigma)} = \frac{1}{1 + \sigma - g_e (x; \sigma)}
\eeq
\cite[\S 4.3]{R62} and \cite{AW94}.  Riordan \cite[(18a), p. 49]{R62} and Tak\'{a}cs \cite{T62} develop a nice recursion for the moments,
 which as a function of $\sigma$ becomes
\beql{w2}
E[W^{n-1} (\sigma)] = \frac{\sigma}{n} \sum_{k=2}^{n} \binom{n}{2} g_k E[W^{n-k} (\sigma)] \qandq E[W^{0} (\sigma)] \equiv 1,
\eeq
where $g_k$ is the $k^{\rm th}$ moment of the service-time cdf.  Clearly, \eqn{w2} can be the source of many integer sequences.

To illustrate,
we now consider the Catalan numbers as service-time moments, which is legitimate by Theorem \ref{thCat}.  We exploit the following result, obtained by
combining \cite[Corollaries 1.3.2 and 1.5.1]{AW87} with Theorem \ref{thCat}.
\begin{lemma} Let $h_1$ be the density of the RBM first moment function, which has the Catalan numbers as its mixing moments.  Then
\beql{h5}
\hat{h}_{1,e} {s} = \hat{h}_1 (s)^2, \quad \mbox{so that} \quad c_e(x) = c(x)^2,
\eeq
where $c(x)$ is the generating function of the Catalan numbers.
\end{lemma}

\begin{theorem}
If the service time pdf is the RBM first moment pdf $h_1$ characterized by its LT in {\em \eqn{h2}}, which has the Catalan numbers as its moments,
then the waiting time mgf in {\em \eqn{w1}} becomes
\beql{rep1}
w(x; \sigma) = \frac{1}{1 + \sigma - \sigma c(x)^2}
\eeq
and the moments of the stationary waiting time pdf are
\beql{rep2}
E[W^{n-1} (\sigma)] = \sigma \sum_{k = 2}^{n} \frac{(n-1)!}{(n-k)!} C_k E[W^{n-k} (\sigma)].
\eeq
Hence, if $\sigma$ is an integer, then $\{E[W^n (\sigma)]\}$ is an integer sequence.
\end{theorem}

For the case $\sigma = 1$, we find $\{E[W^n (1)]\} = 1, 2, 18, 252, 4776, \ldots$, which is not in the OEIS.
Riordan \cite[(19), p. 50]{R62} shows how the moments $E[W^n (\sigma)]$ can be expressed in terms of the multivariate Bell polynomials.


Let $w_n (\sigma) \equiv E[W^n (\sigma)]/n!$ be the associated mixing moments.  From \eqn{rep2}, we obtain a recursion for $w_n (\sigma)$.
\begin{corollary}\label{cor7}
A recursion for the mixing moments $w_n (\sigma)$ defined above is
\beql{rep3}
w_{n-1} (\sigma) = \sigma \sum_{k=2}^{n} C_k w_{n-k} (\sigma), \quad n \ge 1.
\eeq
\end{corollary}
For $\sigma = 1$, we get
 $\{w_n (1)\} = 1, 2, 9, 42, 199, \ldots$,
 which also is not yet in OEIS.

 \subsection{General Birth-and-Death Processes}

 The most relevant previous work connecting queueing theory to integer sequences and providing a suitable framework for generalization
 evidently is the previous work connecting birth-and-death processes to continued fractions; see \cite{AW99a, FG00} and references therein;
 continued fractions are known to be intimately connected to integer sequences.
 We have already mentioned that the stochastic process $Q \equiv \{Q(t): t\ge 0\}$, representing the number of customers in
 an $M/M/1$ queueing model at time $t$, is a birth-and-death stochastic process.  Indeed, it is a special birth-and-death process with
 constant birth rate $\lambda$ and constant death rate $\mu$.  More generally, these rates are functions of the state \cite[\S 6.3]{R07};
 $\lambda_k$ is the rate up and $\mu_k$ is the rate down when $Q(t) = k$.  The general birth-and-death process represents a more general queueing model.

 For a general birth-and-death process $X \equiv \{X(t): t \ge 0\}$, an important quantity is the probability
 \beql{return}
 P_{0,0} (t) \equiv P(X(t) = 0| X(0) = 0), \quad t \ge 0;
 \eeq
 see \cite[\S 11]{AW88a}, \cite[\S 3.1]{FG00} and \cite{K72}.
From \cite[(4.1), p. 771]{FG00}, the LT of $P_{0,0} (t)$, $\hat{P}_{0,0} (s)$, has a
representation as an $S$ fraction
\begin{eqnarray}
\hat{P}_{0,0} (s) & =  &  \frac{1}{s+} \frac{\lambda_0}{1+} \frac{\mu_1} {s+} \frac{\lambda_1}{1+} \frac{\mu_2}{s+} \ldots \nonumber \\
                  & =  &  \frac{z}{1+} \frac{\lambda_0 z}{1+} \frac{\mu_1 z} {1+} \frac{\lambda_1 z}{1+} \frac{\mu_2 z}{1+} \ldots
\end{eqnarray}
where $z \equiv 1/s$.
(The last line follows from \cite[(1.5)]{AW99a} using the sequence $\{c_n: n \ge 0\}$ with $c_{2n} \equiv 1$ and $c_{2n+1} \equiv z$, $n \ge 1$.)
We can find the corresponding power series
\beq
\hat{P}_{0,0} (s) = z(1 - p_1 z + p_2 z^2 - p_3 z^3 + \ldots );
\eeqno
see \cite[(3.9), (3.10), p. 397]{AW99a}.  The sequence $\{1, p_1, p_2, p_3, \ldots\}$ may be an integer sequence.

For the special case of the $M/M/1$ queue, with our scaling we have $\lambda_i = \sigma$ and $\mu_{i+1} = 1 + \sigma$ for all $i$, $i \ge 0$,
so that
\beq
\hat{P}_{0,0} (s) = z \left(\frac{1}{1+} \frac{\sigma z}{1+} \frac{(1 + \sigma) z} {1+} \frac{\sigma z}{1+} \frac{(1 + \sigma) z}{1+} \ldots \right).
\eeqno
Then from Table \ref{tb1} we have the relation
\beql{return2}
\hat{P}_{0,0} (s) = z \hat{p} (z),
\eeq
where $\hat{p} (s)$ is the Laplace transform of the equilibrium time to emptiness.





\section{Proofs}\label{secProofs}

In this concluding section we provide proofs for Theorems \ref{th1} and \ref{th2}.

\subsection{Proof of Theorem \ref{th1}.}

We give two proofs.  The first is direct; the second applies \cite{B09a}.

\paragraph{Direct proof of \eqn{b8}.}


From \cite[pp. 232--233]{T67}, after a change of scale and notation, we have
\beq
b_{n+1} = \left(\frac{1 + \sigma}{n+1}\right) \sum_{k=1}^{n} \binom{n+k}{k}\binom{n-1}{k-1} \sigma^{k-1}.
\eeqno
(This equation is also given by \cite[(21), p. 151]{R68} after we make the identification
that $H_n (z) = p_n (\sigma) = \sigma b_{n+1} (\sigma)/(1 + \sigma)$.)  On the right side,
move $(1+ \sigma)$ inside the sum and identify the coefficients of $\sigma^k$, obtaining
\beq
\left(\frac{1}{n+1}\right) \left(\binom{n+k+1}{k+1}\binom{n-1}{k} + \binom{n+k}{k}\binom{n-1}{k-1}\right)
= \binom{n+k}{n-k}\binom{2k}{k}\left(\frac{1}{k+1}\right).~~~\bsq
\eeqno



\paragraph{Direct proof of \eqn{b8a}.}

We follow the technique of \cite[p. 107]{R62} to establish a three-term recursion.  By the successive differentiation
with respect to $x$ in \eqn{b4}, we find that
\beq
\Psi(x)^2 b''(x; \sigma) = (1+\sigma)(1 + 2 \sigma - x)b' (x; \sigma) + (1+\sigma)^2 b(x; \sigma).
\eeqno
In this equation, make the following substitutions:
\begin{eqnarray}
b(x; \sigma) = \sum_{n=0}^{\infty} b_n (\sigma) x^n, \nonumber \\
b'(x; \sigma) = \sum_{n=0}^{\infty} (n+1)b_{n+1} (\sigma) x^n, \nonumber \\
b''(x; \sigma) = \sum_{n=0}^{\infty} (n+1)(n+2) b_{n+2} (\sigma) x^n. \nonumber
\end{eqnarray}
Then collect and equate the coefficients of $x^{n-1}$ and the result follows.~~~\bsq

\paragraph{Application of \cite{B09a}.}

The idea in this second proof is to directly make connection to the Catalan numbers via their generating function $c$ in \eqn{b7}.
Starting from \eqn{b4}, we apply \eqn{be}, \eqn{b10} and the proof of \eqn{b15} in Theorem \ref{th2} below to conclude that
\beql{k1}
b_e (x; \sigma) = \frac{1}{1-x}c\left(\frac{\sigma x}{(1 - x)^2}\right) \qandq p (x; \sigma) = \frac{1}{1+x}c\left(\frac{(1 +\sigma) x}{(1 + x)^2}\right).
\eeq
We then apply \cite[Proposition 15, p.12]{B09a} (with change of variables $y \equiv \sigma x$ and $y \equiv (1+\sigma) x$, respectively) to immediately deduce the conclusions
\beql{k2}
b_{e, n} (\sigma) = \sum_{k=0}^{n} \binom{n+k}{n-k} C_k \sigma^k \qandq p_{n} (\sigma) = \sum_{k=0}^{n} \binom{n+k}{n-k} (-1)^{n+k} C_k (1 +\sigma)^k
\eeq
given in Corollary \ref{cor2} and Theorem \ref{th2}.  We then can apply \eqn{ex1} and \eqn{be} to get the conclusion in Theorem \ref{th1} for $b$ from $b_e$.
By this line of reasoning, we establish Theorems \ref{th1} and \ref{th2} and Corollaries \ref{cor1} and \ref{cor2} all at once,
and have a new perspective on their relationship.

\subsection{Proof of Theorem \ref{th2}.}

\paragraph{Proof of \eqn{b15}.}

From the representation of $p(x; \sigma)$ in \eqn{b10}, we have
\begin{eqnarray}
p(x; \sigma) & = &\frac{1}{1+\sigma} + \frac{\sigma}{1+\sigma}\left(\frac{2}{1 - x + \Psi(x)}\right) \nonumber \\
& = & \frac{1}{1+\sigma} + \frac{\sigma}{1+\sigma}\left(\frac{1 - x - \Psi(x)}{2 \sigma x}\right) \nonumber \\
& = & \frac{1 + x - \Psi (x)}{2 (1 + \sigma) x} = \frac{2}{1 + x + \Psi(x)}. \nonumber
\end{eqnarray}
On the other hand, by \eqn{b7},
\begin{eqnarray}
\frac{c((1+\sigma)x/(1 +x)^2)}{1+x} & = &\frac{2}{1 + x + (1 + x)\sqrt{1 - 4(1+\sigma) x/(1+x)^2}}  \nonumber \\
& = & \frac{2}{1 + x + \Psi(x)}.~~~\bsq \nonumber
\end{eqnarray}

\paragraph{Proof of \eqn{b15a}.}

Combine \eqn{b8} and \eqn{b10}.~~~\bsq

\paragraph{Proof of \eqn{b15b}.}

Combine \eqn{b8a} and \eqn{b10}.~~~\bsq

\paragraph{Proof of \eqn{b16}.}

From \eqn{b4},
\begin{eqnarray}
b(x; \sigma) & = & \frac{1 + 2 \sigma -x - \Psi(x)}{2 \sigma} = \frac{2(1 + \sigma)}{1 + 2\sigma - x + \Psi(x)} \nonumber \\
&=& \frac{1}{1 - x\left([1+x - \Psi(x)]/[2(1 + \sigma)x]\right)} = \frac{1}{1 - x p(x; \sigma)}.~~~\bsq \nonumber
\end{eqnarray}



\section{Acknowledgment}

The second author was supported by NSF Grant CMMI 0948190.





\pagebreak
\begin{thebibliography}{5}

\bibitem{ACLW95}
J. Abate, G. L. Choudhury, D. M. Lucantoni and W. Whitt, Asymptotic analysis of tail probabilities based on the computation of moments,
{\em Ann. Appl. Prob.} {\bf 5} (1995), 983--1007.



\bibitem{AW87}
J. Abate and W. Whitt, Transient behavior of regulated Brownian motion,
{\em Adv. Appl. Prob.} {\bf 19} (1987), 560--598.


\bibitem{AW88a}
J. Abate and W. Whitt, Transient behavior of the $M/M/1$ queue via Laplace transforms,
{\em Adv. Appl. Prob.} {\bf 20} (1988), 145--178.

\bibitem{AW88b}
J. Abate and W. Whitt, Simple spectral representations for the $M/M/1$ queue,
{\em Queueing Systems} {\bf 3} (1988), 321--346.

\bibitem{AW94}
J. Abate and W. Whitt, Transient behavior of the $M/G/1$ workload process,
{\em Operations Research} {\bf 42} (1994), 750--764.

\bibitem{AW95}
J. Abate and W. Whitt, Limits and approximations for the busy-period distribution in single-server queues,
{\em Prob. Engr. Inform. Sci.} {\bf 9} (1995), 581--602.


\bibitem{AW96}
J. Abate and W. Whitt, An operational calculus for probability distributions via Laplace transforms,
{\em Adv. Appl. Prob.} {\bf 28} (1996), 75--113.


\bibitem{AW99a}
J. Abate and W. Whitt, Computing Laplace transforms for numerical inversion via continued fractions,
{\em INFORMS J. Computing} {\bf 11} (1999), 394--405.

\bibitem{AW99b}
J. Abate and W. Whitt, Explicit $M/G/1$ waiting-time distributions for a class of long-tail service-time distributions..
{\em Operations Research Letters} {\bf 25} (1999), 25--31.

\bibitem{A65}
N. I. Akhiezer, {\em The Classical Moment Problem and Some Related Questions in Analysis}, Oliver and Boyd, London, 1965.

\bibitem{A04}
M. H. Albert, R. E. L. Aldred, M. D. Atkinson, H. P. van Ditmarsch, C. C. Handley
and D. A. Holton, Restricted permutations and queue junping, {\em Discrete Math.} {\bf 287} (2004), 129--133.

\bibitem{AB94}
M. D. Atkinson, R. Beals, Priority queues and permutations,  {\em SIAM J. Comput.} {\bf 23} (1994), 1225--1230.

\bibitem{B09}
P. Barry, {\em A Study of Integer Sequences, Riordan Arrays, Pascal-like Arrays and Hankel Transforms},
PhD dissertation Department of Mathematics, University College Cork, Ireland, 2009.

\bibitem{B09a} P. Barry,  Continued fractions and
transformations of integer sequences,
{\em J. Integer Sequences} {\bf 12} (2009), 
\href{http://www.cs.uwaterloo.ca/journals/JIS/VOL12/Barry3/barry93.html}{Article 09.7.6}.


\bibitem{CBVWJ08}
A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland and W. B. Jones,
{\em Handbook of Continued fractions for Special Functions}, Springer, 2008.

\bibitem{CL96}
G. L. Choudhury and D. M. Lucantoni, Numerical computation of the moments of a probability distribution
from its transform, {\em Oper. Res.} {\bf 44}  (1996), 368--381.

\bibitem{DS97}
H. Dette and W. J. Studden, {\em The Theory of Canonical Moments with Applications to Statistics, Probability and Analysis},
Wiley, 1997.

\bibitem{FG00}
P. Flajolet and F. Guillemin, The formal theory of birth-and-death processes, lattice path combinatorics and continued fractions,
{\em Adv. Appl. Prob.} {\bf 32} (2000), 750--778.

\bibitem{K72}
J. F. C. Kingman, {\em Regenerative Phonomena},
Wiley, 1972.

\bibitem{L01}
J. W. Layman, The Hankel transform and some of its properties,
{\em J. Integer Sequences} {\bf 4} (2001), 
\href{http://www.cs.uwaterloo.ca/journals/JIS/VOL4/LAYMAN/hankel.html}{Article 01.1.5}.


\bibitem{P02}
K. A. Penson, An integral representation of $C_n$,
commentary on (\seqnum{A000108}), in \cite{S10}.

\bibitem{R58}
J. Riordan, {\em An Introduction to Combinatorial Analysis},
 Wiley, 1958.

 \bibitem{R62}
J. Riordan, {\em Stochastic Service Systems},
 Wiley, 1962.

 \bibitem{R68}
J. Riordan, {\em Combinatorial Identities},
 Wiley, 1968.

  \bibitem{R07}
S. M. Ross, {\em Introduction to Probability Models}, ninth edition,
Academic Press, 2007.

 \bibitem{S10}
N. J. A. Sloane,
\href{http://www.research.att.com/~njas/sequences/}{The On-Line
Encylopedia of Integer Sequences} (OEIS), 2010.

 \bibitem{S05}
R. P. Stanley,
Queue problems revisited, Suomen Tehtavaniekat (Proceedings of the Finnish Chess Problem Society) {\bf 59} (2005), 193--203.

 \bibitem{T62}
L. Tak\'{a}cs, {\em Introduction to the Theory of Queues}, Oxford University Press, 1962.

 \bibitem{T67}
L. Tak\'{a}cs,
{\em Combinatorial Methods in the Theory of Stochastic Processes},  Wiley, 1967.

 \bibitem{T90}
 L. Tak\'{a}cs, On the number of distinct forests,  {\em SIAM J. Discrete Math.} {\bf 3} (1990), 574--581.

\bibitem{W85}
W. Whitt,
 The renewal-process stationary-excess operator. {\em J. Appl. Prob.} {\bf 22} (1985), 156--167.

\bibitem{W02}
W. Whitt, {\em Stochastic-Process Limits}, Springer, 2002.


\end{thebibliography}

\bigskip
\hrule
\bigskip

\noindent 2000 {\it Mathematics Subject Classification}:
Primary 11B25; Secondary 60K25.

\noindent \emph{Keywords: }  integer sequence, moment sequence,
probability, probability density function, probability moment
sequence, the moment problem, queue, queueing, Laplace
transforms, moment generating functions, generating functions,
continued fractions, Hankel transform, renewal theory, birth-and-death process,
stationary excess distribution, stationary lifetime
distribution, exponential mixture, single-server queue, $M/M/1$
queue, $M/G/1$ queue, busy period, equilibrium time to
emptiness, waiting time.

\bigskip
\hrule
\bigskip

\noindent  (Concerned with sequences
\seqnum{A000108},
\seqnum{A001003},
\seqnum{A006318},
\seqnum{A093252},
\seqnum{A093526},
\seqnum{A103210},
\seqnum{A103211},
\seqnum{A107841},
\seqnum{A122525},
\seqnum{A131763},
\seqnum{A137216}, and
\seqnum{A155069}.)

\bigskip
\hrule
\bigskip

\vspace*{+.1in} \noindent Received February 19 2010; revised version
received May 2 2010.  
Published in {\it Journal of Integer Sequences}, May 3 2010.

\bigskip
\hrule
\bigskip

\noindent
Return to
\htmladdnormallink{Journal of Integer Sequences home page}{http://www.cs.uwaterloo.ca/journals/JIS/}.
\vskip .1in


\end{document}

.
