††thanks:

# Holonomic Bessel modules and generating functions

Yik Man Chiang<EMAIL_ADDRESS>Department of Mathematics, The Hong Kong
University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong SAR.
,  Avery Ching<EMAIL_ADDRESS>Department of Mathematics, University
of Dundee, Nethergate, Dundee DD1 4HN, Scotland, UK.  and  Xiaoli Lin
<EMAIL_ADDRESS>Department of Mathematics, The Hong Kong University of
Science and Technology, Clear Water Bay, Kowloon, Hong Kong SAR.

###### Abstract.

We have solved a number of holonomic PDEs derived from the Bessel modules
which are related to the generating functions of classical Bessel functions
and the difference Bessel functions recently discovered by Bohner and Cuchta.
This $D$-module approach both unifies and extends generating functions of the
classical and the difference Bessel functions. It shows that the algebraic
structures of the Bessel modules and related modules determine the possible
formats of Bessel’s generating functions studied in this article. As a
consequence of these $D$-modules structures, a number of new recursion
formulae, integral representations and new difference Bessel polynomials have
been discovered. The key ingredients of our argument involve new transmutation
formulae related to the Bessel modules and the construction of $D$-linear maps
between different appropriately constructed submodules. This work can be
viewed as $D$-module approach to Truesdell’s $F$-equation theory specialised
to Bessel functions. The framework presented in this article can be applied to
other special functions.

###### keywords:

$D$-modules, Bessel modules, Generating functions, Difference Bessel
polynomials, Glaisher’s modules, difference Bessel polynomials.

###### 2020 Mathematics Subject Classification:

Primary 32S40, 33E99, 30D10; Secondary 47B37 , 47B47, 12H05, 12H10, 13N10.

###### Contents

  1. 1 Introduction
    1. 1.1 Bessel’s modules
    2. 1.2 Glaisher’s modules
    3. 1.3 Bessel polynomial modules
  2. 2 Preliminaries
    1. 2.1 Linear maps between $D$-modules
    2. 2.2 Integrability criteria
    3. 2.3 Banach algebra and Poincaré-Perron theory
  3. 3 Bessel modules and Bessel functions
    1. 3.1 Weyl series solutions
    2. 3.2 Transmutation formulae
    3. 3.3 Bessel modules
    4. 3.4 Characteristics of Bessel modules
    5. 3.5 Generating function of classical Bessel functions
    6. 3.6 Generating function of difference Bessel functions
    7. 3.7 Integral representations
  4. 4 Half-Bessel modules I: Bessel polynomials
    1. 4.1 Transmutation formulae
    2. 4.2 Reverse Bessel polynomial modules
    3. 4.3 Characteristics of modified half-Bessel modules
    4. 4.4 Bessel polynomials module
    5. 4.5 Characteristic of Bessel polynomial module
  5. 5 Newton transformations
  6. 6 Poisson-type generating functions for (difference) Bessel polynomials
    1. 6.1 Classical reverse Bessel polynomials $\theta_{n}(x)$
    2. 6.2 Difference reverse Bessel polynomials $\theta^{\Delta}_{n}(x)$
    3. 6.3 Classical Bessel polynomials $y_{n}(x)$
    4. 6.4 Difference Bessel polynomials $y^{\Delta}_{n}(x)$
  7. 7 Half-Bessel modules II: Glaisher’s generating functions
    1. 7.1 Characteristics of Glaisher modules
    2. 7.2 Glaisher’s generating functions
    3. 7.3 Difference Glaisher’s generating functions
  8. 8 Discussion
  9. 9 Conclusion
  10. 10 Outlook
  11. A Proofs of some theorems/propositions
    1. A.1 Proof of Proposition 3.33
    2. A.2 Proof of Proposition 4.12
    3. A.3 Proof of Theorem 4.34
    4. A.4 Completion of proof of Theorem 6.11
    5. A.5 Proof of Theorem 6.12
  12. B Difference operators with arbitrary step sizes
  13. C List of tables
    1. C.1 List of commonly used $D$-modules
    2. C.2 List of transmutation formulae
    3. C.3 List of holonomic modules
    4. C.4 List of formulae in $\mathcal{O}_{d}$, $\mathcal{O}_{\Delta}$
    5. C.5 List of PDE systems in complex domains
    6. C.6 List of generating functions

##  1\. Introduction

Generating functions play important roles in both pure and applied mathematics
[27, 67]. In particular, generating functions of Bessel functions play pivotal
roles in the mathematical formulation of electromagnetic wave and acoustic
scattering problems, see for examples, [7, 23]. Other application includes the
study of Rogers-Ramanujan identities [24, 25]. In this paper, we present a
systematic study of generating functions of Bessel functions based on the
Bessel module (1.2) defined below. This is the first of a series of papers
originates from our study of generating functions of classical special
functions using $D$-module methodology. It is known that $D$-modules is an
efficient language in exhibiting algebraic structures and in carrying out
explicit computation of identities of special functions in general, see for
examples [14, 15, 24, 25, 38, 39, 40, 46, 54, 68, 70]. We demonstrate that
this is indeed the case that the $D$-modules approach not only allows us to
better understand some of the classical formulae about special functions but
also to derive their new difference analogues. This work can be considered as
a continuation of the previous works of [43, 44, 51, 60, 58, 61, 62] and to a
less extent of [8, 24, 25] by algebraic methods. A noteworthy feature of this
study is that the main results are first obtained in the $D$-modules framework
before their $D$-modules structures are being manifested in complex domains
differently. Thus the study unifies the classical Bessel functions and the
recently discovered difference Bessel functions amongst other possible
manifestations. The fact that the main results have been derived entirely
within the appropriate $D$-modules that is subject to different manifestations
in the complex plane is a vindication of the viewpoint from umbral calculus
discovered centuries ago.

###  1.1. Bessel’s modules

We extend Bessel’s classical generating function111It was known to Hansen in
1845 and earlier to Jacobi in the 1836 in a modified form [64, §2.1].

$e^{\frac{x}{2}\big{(}t-\frac{1}{t}\big{)}}=\sum_{n=-\infty}^{\infty}J_{n}(x)\,t^{n}$
(1.1)

for the Bessel coefficients $(J_{n}(x))$ which denotes the bilateral sequence
of Bessel functions of the first kind of integer orders, to the case of non-
integer orders $(J_{\nu+n}(x))$ where $\nu\not=0$, amongst other results
obtained, including their difference analogues, by studying the quotient
module

$\mathcal{B}_{\nu}:=\frac{\mathcal{A}_{2}}{\mathcal{A}_{2}(X_{1}\partial_{1}+(\nu+X_{2}\partial_{2})-X_{1}X_{2})+\mathcal{A}_{2}(X_{1}\partial_{1}-(\nu+X_{2}\partial_{2})+{X_{1}}/{X_{2}})},$
(1.2)

be called the Bessel module of order $\nu\in\mathbb{C}$ in this article, where
the $\mathcal{A}_{2}$ denotes the $D$-modules of the Weyl-algebra generated by
$\partial_{i},\,X_{i}\ (j=1,\,2)$ subject to the standard commutation
relations

$[\partial_{j},X_{k}]=\delta_{j,k},\qquad[\partial_{j},\partial_{k}]=0,\qquad[X_{j},X_{k}]=0,\qquad\mbox{for
}j,\,k=1,\,2.$ (1.3)

Here the two relations

$X_{1}\partial_{1}+\nu+X_{2}\partial_{2}-X_{1}X_{2},\qquad
X_{1}\partial_{1}-(\nu+X_{2}\partial_{2})+{X_{1}}/{X_{2}},$ (1.4)

which define the Bessel module $\mathcal{B}_{\nu}$, are abstractions of the
well-known recursions 333Lommel [43, 1868].

$\begin{split}&xJ^{\prime}_{\nu+n}(x)-(\nu+n)J_{\nu+n}(x)=-xJ_{\nu+n+1}(x),\\\
&xJ^{\prime}_{\nu+n}(x)+(\nu+n)J_{\nu+n}(x)=xJ_{\nu+n-1}(x),\end{split}$ (1.5)

about Bessel functions, known as Gauss’ contiguous relations444They and other
generalised ones will be called PDEs in a board sense in this article. The
list of all holonomic PDEs are listed in Appendix §C.3. in special functions,
which hold for every $n\in\mathbb{Z}$, where the $\partial_{1}$ and
$\partial_{2}$ in (1.4) play the roles of differentiation with respect to $x$
and shift of subscripts $n$ respectively. As we shall see that the
$\mathcal{B}_{\nu}$ is holonomic in the sense of Bernstein-Kashiwara-Sato’s
holonomic PDEs theory (see, e.g., [17, 38]), and that solving
$\mathcal{B}_{0}$ gives raise to the (1.1) as a special case of $\nu=0$. The
two relations from (1.4) forms a system of two PDEs from this viewpoint. Thus,
we have derived

$t^{-\nu}e^{\frac{x}{2}\big{(}t-\frac{1}{t}\big{)}}\sim\sum_{n=-\infty}^{\infty}J_{\nu+n}(x)\,t^{n},$
(1.6)

for general $\nu\not=0$, and its difference analogue

$e^{i\pi\nu}\frac{\sin(x-\nu)\pi}{\sin(\pi
x)}\,t^{-\nu}\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}\sim\sum_{n=-\infty}^{\infty}J^{\Delta}_{n+\nu}(x)\,t^{n},$
(1.7)

for the difference Bessel functions $(J_{\nu+n}^{\Delta})$ discovered by
Bohner and Cuchta [4] to be discussed below. Here the “$\sim$” signs used
above denote the Borel-resummation on the right-sides of the corresponding
formulae since both the sequences $(J_{\nu+n})$ and the $(J^{\Delta}_{\nu+n}$)
are $1$-Gervey666See Theorem 3.22 and Theorem 3.39.777For an up-to-date
discussion about the resurgence theory, please see [49, Part II], [59] and
[19].  in the sense that in each compact set, $|J_{\nu+n}(x)|\leq C^{n}n!$ and
$|J^{\Delta}_{\nu+n}(x)|\leq D^{n}n!$ for some positive constants $C,\ D$, see
Example 2.7. Both sides of (1.5) and (1.7) satisfy the systems of PDEs888We
use the terminology “PDEs” in a generalised sense that may include difference
or shift operators or simply members in the $\mathcal{A}_{2}$ in this paper.

$\begin{split}&y_{x}(x,t)+\big{(}1/t-t\big{)}/2\,y(x,t)=0,\\\ &ty_{t}(x,t)+\nu
y(x,t)-{x}/{2}\,\big{(}1/t+t\big{)}\,y(x,t)=0,\end{split}$ (1.8)

and

$\begin{split}&y(x+1,t)-y(x,t)+(1/t-t)/2y(x,t)=0,\\\ &ty_{t}(x,t)+\nu
y(x,t)-x(1/t+t)/2y(x-1,t)=0\end{split}$ (1.9)

respectively.

The method of derivation of the above results and other results in this paper
is based on construction of $D$-linear maps from the Bessel modules
$\mathcal{B}_{\nu}$ (and other $D$-modules) to appropriately manifested
$D$-module analytic function spaces of two variables (e.g.,
$\mathcal{O}_{dd},\ \mathcal{O}_{\Delta d}$). Other quotient $D$-modules
related to the half-Bessel module $\mathcal{B}_{1/2}$ studied in this paper
includes Bessel polynomial modules $\Theta,\mathcal{Y}$ and Glaisher modules
$\mathcal{G}_{\pm}$ which also are listed in Appendix §C.3. This paper can be
considered as a $D$-module study of Truesdell’s $F$-equation theory [62],
which summarised and extended research works about generating functions of
predecessors from the nineteenth century, when applied to Bessel’s generating
functions.999The authors only got to know Truesdell’s work after obtaining the
main results of this paper. However, unlike the analytic approach used in
[62], the existence and uniqueness of the $F$-equation theory are replaced by
the holonomicity of the $D$-modules concerned. The paper aims to show that
$D$-modules is not only a natural language but also an efficient computational
tool to study special functions. We have built up our applications of
$D$-modules that represent some of the simplest holonomic equations which give
raise to the exponential, trigonometric functions, etc from scratch in order
to deal with the quotient modules studied in this paper. This is because we
could reach our main results efficiently with our self-contained approach
instead of utilising the full force of $D$-modules theory from, for examples,
[38, 56]. This is especially the case when we derive difference analogues of
the classical results. Closer to our approach but differ substantially in
details were the pioneering works of Wilf and Zeilberger [70, 68] in which the
authors used holonomic $D$-modules methodology to derive hypergeometric type
identities. Our approach is also different from the Lie algebraic approach of
generating functions used by Weisner [65], see also [47]. The choice to start
with Bessel functions is a natural one in regard to them being ubiquity in
many research areas and wide range of applications, see for examples [2, 26,
43, 63, 66, 64].

For easy of computation of the relevant generating functions, we add and
subtract the two relations (1.4) to obtain, respectively, the two alternative
relations

$2X_{1}\partial_{1}-X_{1}(X_{2}-1/X_{2}),\qquad
2\nu+2X_{2}\partial_{2}-X_{1}(X_{2}+1/X_{2})$ (1.10)

for which we can, without loss of generality, eliminate the term “$X_{1}$”
altogether from the first relation above. As an example with $\nu=0$, after a
“change of variables” 101010This change of variables corresponds to the
“characteristic method”., the (1.10) is rewritten as

$2\hat{X}_{1}\hat{\partial}_{1}-\hat{X}_{1},\quad
2\hat{X}_{2}\hat{\partial}_{2}+\hat{X}_{2},$

where the
$\\{\hat{\partial}_{1},\hat{X}_{1},\hat{\partial}_{2},\hat{X}_{2}\\}$ satisfy
similar commutation relations as those in the (1.3), and which when solved
simultaneously leads us to the solution

${\mathop{\rm E}}\Big{[}\frac{X_{1}}{2}(X_{2}-\frac{1}{X_{2}})\Big{]},$ (1.11)

where the ${\mathop{\rm E}}$ denotes the Weyl exponential (see Example 3.3).
Right multiplication by $e^{\frac{x}{2}(x-1/x)}$ gives the map
$\mathcal{B}_{0}\to\mathcal{O}_{dd}$ where the $\mathcal{O}_{dd}$ has the
manifestation, as it is called in this paper,

$\begin{aligned} &(\partial_{1}f)(x,t)=f_{x}(x,t),\\\
&(X_{1}f)(x,t)=xf(x,t),\end{aligned}\qquad\begin{aligned}
&(\partial_{2}f)(x,t)=f_{t}(x,t),\\\ &(X_{2}f)(x,t)=tf(x,t)\end{aligned}$
(1.12)

so that the $\mathcal{O}_{dd}$ becomes an $\mathcal{A}_{2}$-module, and the
map itself is clearly well-defined. It turns out that the system of PDEs
(1.10) become, when $\nu=0$, the following system of PDEs

$\begin{split}&y_{x}(x,t)+\big{(}1/t-t\big{)}/2\,y(x,t)=0,\\\
&ty_{t}(x,t)-x\,\big{(}1/t+t\big{)}/2\,y(x,t)=0,\end{split}$ (1.13)

via an $A_{2}$-linear map, and that both sides of the (1.1) satisfy. Indeed,
the $e^{\frac{x}{2}\big{(}t-\frac{1}{t}\big{)}}$ is an image of (1.11) via
$\mathcal{A}_{2}/(\mathcal{A}_{2}\partial_{1}+\mathcal{A}_{2}\partial_{2})\to\mathcal{O}_{dd}$.
On the other hand, we may equip the $\mathcal{O}^{\mathbb{Z}}$ consisting of
sequences of analytic functions $(f_{n})_{-\infty}^{\infty}$ all defined on an
appropriate domain with a $D$-module structure by the manifestation

$\begin{aligned} &(\partial_{1}f)_{n}(x)=f^{\prime}_{n}(x),\\\
&(X_{1}f)_{n}(x)=xf_{n}(x),\end{aligned}\qquad\begin{aligned}
&(\partial_{2}f)_{n}(x)=(n+1)f_{n+1}(x),\\\
&(X_{2}f)_{n}(x)=f_{n-1}(x).\end{aligned}$ (1.14)

It is easily verified that the $\mathcal{O}^{\mathbb{Z}}$ becomes an
$\mathcal{A}_{2}$-module which is denoted by $\mathcal{O}^{\mathbb{Z}}_{d}$,
so that right multiplication by $(J_{n})$ shows that the map
$\mathfrak{j}:\mathcal{B}_{0}\to\mathcal{O}^{\mathbb{Z}}_{d}$ is also
$\mathcal{A}_{2}$-linear. Moreover, we have the systems111111According to
Watson [64, p. 18, p. 45] these formulae were discovered by Bessel in 1824 and
later the (1.5)by Lommel in 1868 [43] for general $\nu\not=0$.

$\begin{split}&xJ^{\prime}_{n}(x)-nJ_{n}(x)=-xJ_{n+1}(x),\\\
&xJ^{\prime}_{n}(x)+nJ_{n}(x)=xJ_{n-1}(x),\end{split}$ (1.15)

that lie in the kernel of the $\mathcal{A}_{2}$-linear map
$\mathfrak{j}:\mathcal{B}_{0}\to\mathcal{O}^{\mathbb{Z}}_{d}$, which are
special cases of the (1.5) when $\nu=0$. Then following diagram commutes

${\mathcal{B}_{0}}$${\mathcal{O}^{\mathbb{Z}}_{d}}$${\widetilde{\mathcal{A}}_{2}}$${\mathcal{O}_{dd},}$$\scriptstyle{\mathfrak{j}}$$\scriptstyle{\times{\mathop{\rm
E}}[\frac{X_{1}}{2}(X_{2}-\frac{1}{X_{2}})]}$$\scriptstyle{\mathfrak{z}}$$\scriptstyle{\times
1}$ (1.16)

where
$\widetilde{\mathcal{A}}_{2}:=\overline{\mathcal{A}_{2}/\big{[}\mathcal{A}_{2}\partial_{1}+\mathcal{A}_{2}\partial_{2}\big{]}}$
is the completion of
$\mathcal{A}_{2}/\big{[}\mathcal{A}_{2}\partial_{1}+\mathcal{A}_{2}\partial_{2}]$
and where the map $\mathfrak{z}$ is the $z$-transform, from which the
classical formula (1.1) essentially follows because of the holonomicity of
$\mathcal{B}_{0}$. A more general consideration gives the formula (1.6) for
$\mathcal{B}_{\nu}$ (Theorem 3.24).

The above description illustrates how the $D$-module approach allows us to see
each generating function is a solution for the different holonomic $D$-modules
studied (see the table in Appendix C.3). They are obtained by solving the
systems of holonomic PDEs from, for examples, the (1.9), (1.13), (1.22),
(1.25), (1.27) which are different manifestations inherited from the same
holonomic modules listed in Appendix C.3.

Before we discuss the system (1.9), let us introduce the difference Bessel
function of order $\nu$ recently discovered by Bohner and Cuthta [4]

$J^{\Delta}_{\nu}(x):=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{2^{\nu+2k}k!\Gamma(\nu+k+1)}(x)_{\nu+2k},$
(1.17)

where $(x)_{\nu}=\Gamma(x+1)/\Gamma(x+1-\nu)$ for an arbitrary $\nu$ from our
$D$-module perspective. Bohner and Cuthta [4] shows that the Newton series
solves the difference Bessel equation

$x(x-1)\triangle^{2}y(x-2)+x\triangle y(x-1)+x(x-1)y(x-2)-\nu^{2}y(x)=0.$
(1.18)

The first two authors and Tsang have shown in [13] that the function (1.17)
and the equation (1.18) are manifestations of the Weyl-algebraic Bessel (3.22)
and the Weyl-Bessel operator

$(X\partial)^{2}+(X^{2}-\nu^{2})$

in the analytic function space $\mathcal{O}_{\Delta}$, where the subscript
$\Delta$ indicates that the $\partial$ in $\mathbb{C}\langle\partial,X\rangle$
subjects to $[\partial,X]=\partial X-X\partial=1$, that is in
$\mathcal{A}_{1}$, is interpreted as a forward difference operator instead,
that is, $\partial f(x)=\Delta f(x)=f(x+1)-f(x)$ and $Xf(x)=xf(x-1)$ (see
Example 3.32.). Indeed, one recovers the classical Bessel equation

$x^{2}y^{\prime\prime}(x)+xy^{\prime}(x)+(x^{2}-\nu^{2})y(x)=0.$ (1.19)

and classical Bessel function with the well-known manifestation $\partial
f(x)=f^{\prime}(x),\ Xf(x)=xf(x)$.

Historically, both Sonine [60] and Nielsen [51] (see also [64, §3.9]) amongst
others from the second half of nineteenth century studied Bessel functions
from the system of partial differential equations (1.5) instead of (1.19).
Indeed, Nielson distinguished the solutions to the system (1.5) called
cylindrical functions against those from the (1.19). These earlier studies
were in line with Truesdell’s “$F$-equation theory” [62] as noted by Truesdell
himself, that we have already mentioned.

An advantage of our $D$-modules approach to the generating function (1.6) is
that we could modify the $\mathcal{A}_{2}$-module structures (1.12) and
(LABEL:E:manifestation_2), which represent one such manifestation that we have
mentioned above, into different ones, say to

$(\partial_{1}f)(x,t)=\Delta_{x}f(x,t)=f(x+1,t)-f(x,t),\quad(X_{1}f)(x,t)=xf(x-1,t)$
(1.20)

and

$(\partial_{1}f)_{n}(x,t)=\Delta_{x}f_{n}(x,t)=f_{n}(x+1,t)-f_{n}(x,t),\quad(X_{1}f)_{n}(x,t)=xf_{n}(x-1,t)$

on $\mathcal{O}_{\Delta d}$ and $\mathcal{O}^{\mathbb{Z}}_{\Delta}$
respectively, while keeping the defining properties of the corresponding
$\partial_{2},\,X_{2}$ in (1.12) and (LABEL:E:manifestation_2) respectively
unchanged. As a result, we arrive at the commutative diagram

${\mathcal{B}_{0}}$${\mathcal{O}^{\mathbb{Z}}_{\Delta}}$${\widetilde{A}_{2}}$${\mathcal{O}_{\Delta
d}}$$\scriptstyle{\mathfrak{j}_{\Delta}}$$\scriptstyle{\times{\mathop{\rm
E}}\big{[}\frac{X_{1}}{2}\big{(}X_{2}-\frac{1}{X_{2}}\big{)}\big{]}}$$\scriptstyle{\mathfrak{z}_{\Delta}}$$\scriptstyle{\times
1}$

similar to the (1.16) for this new $\mathcal{A}_{2}$-module manifestation,
which implies the new generating function

$\Big{[}\frac{1}{2}(t-\frac{1}{t})+1\Big{]}^{x}=\sum_{n=-\infty}^{\infty}J^{\vartriangle}_{n}(x)\,t^{n}$
(1.21)

when $\nu=0$. The full statements will be given in Theorem 3.43 and Theorem
3.38 respectively. Here, both sides of (1.7) satisfy the system of PDEs

$\begin{split}&y(x+1,t)-y(x,t)+(1/t-t)/2y(x,t)=0,\\\
&ty_{t}(x,t)-x(1/t+t)/2y(x-1,t)=0\end{split}$ (1.22)

as images of the (1.10), with $\nu=0$, which are the difference analogues of
those in (1.13) in our $D$-modules interpretation. They are the special cases
of (1.9) when $\nu=0$. The $(J^{\triangle}_{n}(x))$ in (1.21) are precisely
the difference Bessel functions $(J^{\triangle}_{\nu+n}(x))$ discovered by
Bohner and Cuchta [4]. The reason that we can derive the infinite sum
representation of the generating function on the right-hand side of (1.21) is
because of the following formulae

$\displaystyle x\Delta
J^{\Delta}_{n}(x-1)+nJ^{\Delta}_{n}(x)-xJ^{\Delta}_{n-1}(x-1)$
$\displaystyle=0,$ $\displaystyle x\Delta J^{\Delta}_{n}(x-1)-\nu
J^{\Delta}_{n}(x)+xJ^{\Delta}_{n+1}(x-1)$ $\displaystyle=0,$

derived by brute-force verification in [4] are manifestation of the (1.4) in
the kernel of the $\mathcal{A}_{2}$-linear map
$\mathcal{B}_{0}\to\mathcal{O}^{\mathbb{Z}}_{\Delta}$ as similar to the case
of the (1.5).

If we modify the the above commutator to arbitrary step size $h$, i.e.,
$\partial f(x)=\triangle_{h}f(x)=\frac{1}{h}\big{(}f(x+h)-f(x)\big{)}$ and
$Xf(x)=xf(x-h)$, then the generating function (1.21) is replaced by

$\Big{[}\frac{h}{2}(t-\frac{1}{t})+1\Big{]}^{x/h}=\sum_{n=-\infty}^{\infty}J^{\vartriangle_{h}}_{n}(x)\,t^{n}.$
(1.23)

Clearly, we recover the classical formula (1.1) from the (1.23) after letting
$h\to 0$ on both sides of the above formula. See Theorem B.5 in the Appendix A
for more details. By a similar consideration, the formulae (1.5) become

$\begin{split}&x\triangle\mathscr{C}^{\Delta}_{\nu+n}(x-1)+(\nu+n)\mathscr{C}^{\Delta}_{\nu+n}(x)-x\mathscr{C}^{\Delta}_{\nu+n-1}(x-1)=0,\\\
&x\triangle\mathscr{C}^{\Delta}_{\nu+n}(x-1)-(\nu+n)\mathscr{C}^{\Delta}_{\nu+n}(x)+x\mathscr{C}^{\Delta}_{\nu+n+1}(x-1)=0.\\\
\end{split}$ (1.24)

we the difference manifestation of $\mathcal{A}_{2}$. If we choose
$\mathscr{C}_{\nu}(x)=J^{\Delta}_{\nu}(x)$ in the above two formulae for
$J^{\Delta}_{\nu}(x)$, then we recover the two corresponding formulae that
Bohner and Cuthta derived in [4] . Combining the PDEs (1.9) and (1.24) and a
commutative diagram similar to the last commutative diagram above implies the
new generating function (1.7).

The difference Bessel functions are not discussed in Nikifornov et al. We
mention that the difference Bessel functions have not been discussed in recent
major works on difference special functions such as [52] and [53], which
mostly focus on orthogonal polynomials.

Before we discuss similar formulae that we have obtained for the new
difference analogues of Bessel polynomials and Glaisher’s trigonometric type
generating functions below, let us first return to the discussion mentioned in
the first paragraph about the aspect of our Weyl-algebraic methodology in
solving121212See Definition 3.1. the two “PDEs” (1.4) in the context of the
Bessel module $\mathcal{B}_{\nu}$. In fact, the process of a priori
computation for the “generating functions” greatly simplifies the entire
process of finding the generating function that are classically written within
suitably defined analytic function spaces. This a priori computation is
therefore no less significant than the conventional approach done in a complex
domain.

As a by-product of (1.7) of “extracting the residue” from the formula (1.7),
we obtain the integral representation

$J^{\Delta}_{\nu}(x)=e^{i\pi\nu}\dfrac{\sin(x-\nu)\pi}{\sin(\pi
x)}\dfrac{1}{2\pi
i}\int_{\infty}^{(0+)}t^{-\nu-1}\Big{[}\frac{1}{2}(t-\frac{1}{t})+1\Big{]}^{x}dt,\quad\mathop{\rm
Re}(x)<\mathop{\rm Re}(\nu)$

where the integration path is a Hankel-type contour (Theorem 3.45). This
integral representation is the difference analogue of the classical Schläfli-
Sonine integral of Corollary 3.29 [64, p. 176]. Thus, the classical Schläfli-
Sonine integral can be considered as “residue extraction” from the formula
(1.6).

The rest of this paper will focus on two main variants of the Bessel module of
order half $\mathcal{B}_{1/2}$. They are related to classical Glaisher’s
trigonometric type generating functions and the classical Bessel polynomials.
As a result we have also obtained their difference analogues. The solving of
these variants of $\mathcal{B}_{1/2}$ (Theorem 4.22, Theorem 6.8, Theorem 7.4,
and Theorem 7.7) are more complicated in the sense that certain quadratic
extensions on the base Weyl-algebra are needed.

###  1.2. Glaisher’s modules

The treatment of Glaisher’s Poisson type generating functions [64, p. 140]

$\sqrt{\frac{2}{\pi
x}}\cos\sqrt{x^{2}-2xt}=\sum_{n=0}^{\infty}J_{n-\frac{1}{2}}(x)\,\frac{t^{n}}{n!},$

which is valid for $2|t|<|x|$, requires a change of characteristics variables
to (1.4) with $\nu=-\frac{1}{2}$ different from that for obtaining the (1.11),
so that the alternative $\mathcal{A}_{2}$-module, called the negative
Glaisher’s module is given by

$\mathcal{G}_{-}=\dfrac{\mathcal{A}_{2}}{\mathcal{A}_{2}\big{(}W_{1}\Xi_{1}+(W_{2}\Xi_{2}-\frac{1}{2})-W_{1}\Xi_{2}^{-1}\big{)}+\mathcal{A}_{2}\big{(}W_{1}\Xi_{1}-(W_{2}\Xi_{2}-\frac{1}{2})+W_{1}\Xi_{2}\big{)}},$

in which the $\Xi_{1},W_{1}$ and $\Xi_{2},\,W_{2}$ play the roles of
$\partial_{1},\,X_{1}$ and $\partial_{2},\,X_{2}$ respectively, as defined in
(1.3).

As an example, after a “further change of variables” and solving the resulting
PDEs finally results in the following difference analogue of Glaisher’s
generating functions for difference Bessel functions with half-integer orders
(Theorem 7.12):

$\sqrt{\frac{2}{\pi}}\dfrac{e^{-i\pi x}}{2i\sin(\pi
x)}\int_{-\infty}^{(0+)}\frac{e^{\lambda}(-\lambda)^{-x-1}\lambda^{-\frac{1}{2}}\cos\sqrt{\lambda^{2}-2\lambda
t}}{\Gamma(-x)}\,d\lambda=\sum_{n=0}^{\infty}J^{\Delta}_{n-\frac{1}{2}}(x)\frac{t^{n}}{n!},$

valid on the whole $\mathbb{C}\times\mathbb{C}$. In particular, the generating
function is a solution to the system of delay-differential
equations131313These equations are generally called PDEs in a generalised
sense in this paper. (Theorem 7.10)

$\begin{split}&tf_{tt}(x,\,t)+\big{(}x+1/2\big{)}f_{t}(x,\,t)-xf_{t}(x-1,\,t)-xf(x-1,\,t)=0,\\\
&tf_{t}(x,\,t)-xf_{t}(x-1,\,t)+xf(x-1,\,t)-\big{(}x+1/2\big{)}f(x,\,t)=0.\end{split}$
(1.25)

###  1.3. Bessel polynomial modules

Next we consider half Bessel-module $\mathcal{B}_{1/2}$ which is more natural
to treat the classical reverse polynomials $\theta_{n}(x)$ instead of the
$y_{n}(x)=x^{n}\theta_{n}(1/x)$. If $X\partial$ is replaced by
$X\partial-X-1/2$ in (1.4), then the $\mathcal{B}_{1/2}$ is replaced by the
reverse Bessel polynomial module

$\Theta=\dfrac{\mathcal{A}_{2}}{\mathcal{A}_{2}[\partial_{1}-1+X_{1}/\partial_{2}]+\mathcal{A}_{2}[X_{1}\partial_{1}-2X_{2}\partial_{2}-1-X_{1}+\partial_{2}]}.$

such that the two generators $\partial_{1}-1+X_{1}/\partial_{2}$ and
$X_{1}\partial_{1}-2X_{2}\partial_{2}-1-X_{1}+\partial_{2}$ would generate a
third element
$X_{1}\partial_{1}^{2}-2(X_{1}+X_{2}\partial_{2})\partial_{1}+2X_{2}\partial_{2}$
which is a Weyl-algebraic analogue of the classical Bessel polynomial
equation. In particular, we have derived the difference equation (Theorem
4.17)

$(x-2n)\,\theta_{n}^{\Delta}(x+1)-4(x-n)\,\theta_{n}^{\Delta}(x)+3x\,\theta_{n}^{\Delta}(x-1)=0,$
(1.26)

which is satisfied by difference reverse Bessel (Newton) polynomial of degree
$n$

$\theta^{\Delta}_{n}(x):=\sum_{k=0}^{n}\frac{(n+k)!}{2^{k}(n-k)!\,k!}\,(x)_{n-k}.$

In a similar fashion, the two generators $\partial_{1}-1+X_{1}/\partial_{2}$
and $X_{1}\partial_{1}-2X_{2}\partial_{2}-1-X_{1}+\partial_{2}$ implies two
delay-difference equations for difference reversed Bessel polynomials for each
$n$ (Theorem 4.18):

$\theta^{\Delta}_{n}(x+1)-2\theta^{\Delta}_{n}(x)+x\theta^{\Delta}_{n-1}(x-1)=0,$

and

$\theta^{\Delta}_{n+1}(x)+(x-2n-1)\theta^{\Delta}_{n}(x)-2x\,\theta^{\Delta}_{n}(x-1)=0.$

A “characteristic change of variables” to the module $\Theta$ happens to be
$\rho^{2}=1-2X_{2}$ and we formulate a “quadratic extension of $\Theta$” in
Theorem 6.6

$\Theta(\rho):=\dfrac{\mathcal{A}_{2}(\rho)}{\mathcal{A}_{2}(\rho)[-\rho^{2}\partial_{2}^{2}+3\partial_{2}+X_{1}^{2}]+\mathcal{A}_{2}(\rho)[X_{1}\partial_{1}+\rho^{2}\partial_{2}-1-X_{1}]}.$

so that the commutative diagram of $\mathcal{A}_{2}$-linear maps enable us to
derive a generating function for the difference reverse Bessel polynomials

$\frac{e^{-i\pi x}}{2\pi i\sin\pi
x\Gamma(-x)}\int_{-\infty}^{(0+)}e^{\lambda}(-\lambda)^{-x-1}\frac{1}{\sqrt{1-2t}}\exp\big{[}\lambda(1-\sqrt{1-2t})\big{]}\,d\lambda=\sum_{n=0}^{\infty}\theta_{n}^{\Delta}(x)\,\frac{t^{n}}{n!}$

that is valid in $\mathbb{C}\times B(0,\,\frac{1}{2})$ in Theorem 6.8. The
generating function for the difference reverse Bessel polynomials satisfies
the system of delay-differential equations (Theorem 6.15)

$\begin{split}&f_{t}(x+1,\,t)-2f_{t}(x,\,t)+xf(x-1,\,t)=0,\\\
&(1-2t)f_{t}(x,t)+(x-1)f(x,t)-2xf(x-1,t)=0.\end{split}$ (1.27)

This paper is organised as follows. Section 2 explains some examples of
$D$-modules structures on certain spaces of analytic functions and spaces of
sequences of analytic functions $\mathcal{O}^{\mathbb{Z}}$, as well as
$D$-linear maps between them. A list summarising these entities is given in
Appendix C.3. The remaining parts of the §2 gives a rudiment of holonomic
$D$-modules theory needed in this article. A kind of generalised Borel
resummation is also introduced. §3 introduces the Bessel modules
$\mathcal{B}_{\nu}$ for arbitrary $\nu$ and the method of characteristic in
solving a system of PDEs. Before discussing about the Bessel modules, we also
discuss analogues of “elementary functions” such as the Weyl-exponentials,
Weyl-trigonometry, and Weyl-Bessel. This section also contains the statements
and proofs of generating functions for the classical Bessel functions and
difference Bessel functions, plus some immediate applications of these results
to yield difference analogue of the classical Schläfli-Sonine integral
representation from classical Bessel functions to difference Bessel functions.
In section §4 we modify the Bessel module $\mathcal{B}_{1/2}$ to the “reverse
Bessel module” $\Theta$ by the introducing an integrating factor. We work out
a characteristic for $\Theta$ as well as a variant $\mathcal{Y}$ corresponding
to the usual Bessel polynomials $y_{n}$. These characteristics will be applied
to find (Poisson type) generating functions in due course. As applications of
the $D$-modules approach, we derive, in a uniform way, the well-known three-
term recursion formulae for Bessel polynomials as well as the corresponding
ones for the difference reverse Bessel polynomials. Section §5 introduces the
Newton transform, a kind of generalised Mellin transform but with very
different purpose as a $D$-linear map from $\mathcal{O}_{d}$ to
$\mathcal{O}_{\Delta}$, namely the $D$-modules of analytic functions endowed
with a forward difference operator. The Newton transformation serves as an
effective tool to obtain generating functions for difference Bessel functions.
Utilising results obtained from §4, we give detailed proofs of generating
functions for the reverse Bessel polynomials as well as their difference
analogues in §6.2. Section §7 take the half-Bessel modules from §4 a step
further to study trigonometric-type generating functions for the Bessel
functions first obtained by Glaisher [64]. Before the concluding section §9,
§8 discusses our main results against previous results of Truesdell amongst
others with some historical notes. In particular, it focus the philosophy that
Truesdell pursued, and similarly of Wilf and Zeilberger, to find general
methods of discovery of what kinds of formulae should exist rather than by ad
hoc manipulations of formulae, a viewpoint shared by this study. Appendix A
contains proofs of various Theorems and Propositions stated but not fully
verified in all aspects in the main text. Appendix B extends the study of the
manifestation of the Bessel modules $\mathcal{B}_{\nu}$ to
$\mathcal{O}_{\Delta_{h}}$, that is finite difference with arbitrary step size
$h$. In the case of $\mathcal{B}_{0}$, the generating function derived unifies
the case when $h=1$ for the difference Bessel functions and the classical
Bessel functions when $h\to 0$. This presents a further unification of the
generating functions of the classical Bessel functions and that of difference
Bessel functions of arbitrary step sizes. Appendix §C gives lists of
$D$-modules used, transmutation formulae, the systems of PDEs encountered for
various generating functions defined in their appropriate analytic functions
space of two variables.

##  2\. Preliminaries

###  2.1. Linear maps between $D$-modules

######  Definition 2.1 ([13]).

Let $n\in\mathbb{N}.$ The Weyl-algebra $\mathcal{A}_{n}$ is the
$\mathbb{C}$-algebra with $2n$ generators
$X_{1},\ldots,X_{n},\partial_{1},\ldots,\partial_{n}$ subject to the relations

$[X_{i},X_{j}]=0;\quad[\partial_{i},\partial_{j}]=0;\quad[\partial_{i},X_{j}]=\delta_{ij}.$

In this paper, we denote $\mathcal{A}_{1}=\mathbb{C}\langle X,\partial\rangle$
or by $\mathcal{A}$ when this is clear from the context.

We gather here preliminary material about $D$-module that we shall employ to
exhibit a unified theory for the classical Bessel functions and generating
functions for the recently discovered difference Bessel functions by Bohner
and Cuchta [4] in the next sections. More specifically, amongst the generating
functions studied include

  * •

extensions of generating functions of various kinds for the classical Bessel
functions of arbitrary order $\nu$;

  * •

new generating functions for difference Bessel functions [4];

  * •

new generating functions of Glaisher (Poisson)-type for difference Bessel
functions;

  * •

new generating functions for new difference Bessel polynomials.

The primary concern of this paper is on ${D}$-linear maps defined on the
Bessel modules and their induced generating functions to be introduced in the
following two sections. Hence some generating functions found for difference
Bessel functions may not converge in the same way as those of their classical
counterparts.

The principle of obtaining these generating functions behind is based on the
finite dimensionality of the space of $D$-linear maps from $\mathcal{B}_{\nu}$
to different $\mathcal{A}_{2}$-modules. The first step is the introduction of
some examples of $\mathcal{A}_{n}$-modules and linear maps between them which
are frequently used in this article.

######  Definition 2.2.

We denote by $\mathcal{A}(1/X)$ the Weyl-algebra
$\mathbb{C}\langle\partial,X\rangle(1/X)=\mathbb{C}\langle\partial,X,{1}/{X}\rangle$
modulo the aforementioned relations.

######  Example 2.3.

Let $\mathcal{O}$ be the space of analytic functions on a certain domain. Then
the $\mathcal{O}$ becomes a left $\mathcal{A}$-module denoted by
$\mathcal{O}_{d}$, if one endows it with the structure

$\begin{array}[]{l}(\partial f)(t)=f^{\prime}(t),\\\
(Xf)(t)=tf(t)\quad\mbox{for all }f\in\mathcal{O}\mbox{ and }t\mbox{ in this
domain}.\end{array}$

Similarly, let $\mathcal{O}_{2}$ be the space of analytic functions on a
certain domain in $\mathbb{C}\times\mathbb{C}$. Then the $\mathcal{O}_{2}$
becomes a left $\mathcal{A}_{2}$-module denoted by $\mathcal{O}_{dd}$, if one
endows it with the structure

$\begin{split}(\partial_{1}f)(x,\,t)&=f_{x}(x,\,t),\quad(X_{1}f)(x,\,t)=xf(x,\,t),\\\
(\partial_{2}f)(x,\,t)&=f_{t}(x,\,t),\quad(X_{2}f)(x,\,t)=tf(x,\,t).\end{split}$
(2.1)

######  Example 2.4.

Let $\mathcal{O}$ be the space of analytic functions. Then the $\mathcal{O}$
becomes a left $\mathcal{A}$-module, denoted by $\mathcal{O}_{\Delta}$, if one
endows it with the structure

$\begin{split}&(\partial f)(t)=f(t+1)-f(t),\\\ &(Xf)(t)=tf(t-1).\end{split}$

######  Example 2.5.

Let $\mathbb{C}^{\mathbb{Z}}$ be the space of all functions
$\mathbb{Z}\to\mathbb{C}$ (bilateral sequences). Let

$a=(a_{n})=(\cdots,a_{-2},\,a_{-1},\,a_{0},\,a_{1},\,a_{2},\,\cdots)$

denote a bilateral sequence where $a_{k}\in\mathbb{C}$. Then the
$\mathbb{C}^{\mathbb{Z}}$ becomes a left $\mathcal{A}$-module by defining

$\begin{array}[]{l}(\partial a)_{n}=(n+1)a_{n+1}\\\
(Xa)_{n}=a_{n-1}\quad\mbox{ for all bilateral sequences }(a_{n}).\end{array}$

Various growth conditions can be imposed to such bilateral sequences.

####  2.1.1. $z$-transform I

######  Theorem 2.6 ((Generating function of numeric sequences)).

Suppose that $\mathbb{C}^{\mathbb{Z}}$ consists only of sequences
$\\{a_{n}\\}$ satisfying that

$0\leq\limsup_{n\to\infty}{|a_{-n}|^{1/n}}<\frac{1}{\limsup_{n\to\infty}{|a_{n}|^{1/n}}}\leq+\infty,$

and let $\mathcal{O}$ be the space of analytic functions defined on a certain
annulus centred at $0$. Then the z-transform

$\begin{array}[]{rcl}\mathfrak{z}:\mathbb{C}^{\mathbb{Z}}&\longrightarrow&\mathcal{O}_{d},\\\
(a_{n})&\longmapsto&\displaystyle\sum_{n=-\infty}^{\infty}a_{n}t^{n}.\end{array}$

is left $\mathcal{A}$-linear.

###### Proof.

The fact that the map $\mathfrak{z}$ is $\mathbb{C}$-linear is obvious. The
verification that both

$\mathfrak{z}\,\partial(a_{n})=\partial\,\mathfrak{z}(a_{n}),\quad\mbox{and}\quad\mathfrak{z}\,X(a_{n})=X\,\mathfrak{z}(a_{n})$

hold are routine, with the understanding that the maps $\partial$ and $X$ are
interpreted appropriately in $\mathbb{C}^{\mathbb{Z}}$ and $\mathcal{O}_{d}$
respectively. Hence the map $\mathfrak{z}$ is $\mathcal{A}$-linear. Clearly,
it is both injective and surjective. ∎

In fact, Theorem 2.6 can be extended to Theorem 2.6′ below, which handles a
larger class of sequences $\\{a_{n}\\}$ so that
$\sum_{n=-\infty}^{\infty}a_{n}t^{n}$ may not converge.

######  Definition 2.7 ((Borel resummation)).

Suppose that $\\{a_{n}\\}\in\mathbb{C}^{\mathbb{Z}}$ satisfies both

$\limsup_{n\to\infty}{|a_{n}|^{1/n}}<+\infty\quad\mbox{and}\quad\limsup_{n\to\infty}{\left|\frac{a_{-n}}{n!}\right|^{1/n}}<+\infty.$
(2.2)

Then the Borel resummation of the Laurent series
$\sum_{n=-\infty}^{\infty}{a_{n}t^{n}}$ is defined by taking the Laplace
transform of the Borel transform of its principal part, i.e.,

$\int_{0}^{+\infty}{e^{-st}\left(\sum_{n=0}^{\infty}{\frac{a_{-n-1}}{n!}}s^{n}\right)\,ds}+\sum_{n=0}^{\infty}{a_{n}t^{n}}.$
(2.3)

######  Remark 2.8.

A power series or Laurent series whose sequence of coefficients $\\{a_{n}\\}$
satisfy that

$\limsup_{n\to\infty}{\left|\frac{a_{n}}{n!}\right|^{1/n}}<+\infty$

is said to be $1$-Gevrey. Note that the Borel resummation of a $1$-Gevrey
Laurent series always exists as an analytic function, and the Borel
resummation of a convergent Laurent series is just its ordinary sum, which is
again an analytic function. We refer the reader to [49, Part II], [59] and
[19] for general discussions of $1$-Gevrey or $m$-Gevrey series and their
Borel-resummation.

####  2.1.2. $z$-transform II

######  Theorem 2.6′ ((Generating function of numeric sequences)).

Suppose that $\mathbb{C}^{\mathbb{Z}}$ consists only of sequences
$\\{a_{n}\\}$ satisfying both

$\limsup_{n\to\infty}{|a_{n}|^{1/n}}<+\infty\quad\mbox{and}\quad\limsup_{n\to\infty}{\left|\frac{a_{-n}}{n!}\right|^{1/n}}<+\infty,$

and let $\mathcal{O}$ be the space of analytic functions defined on a certain
region. Then the z-transform

$\begin{array}[]{rcl}\mathfrak{z}:\mathbb{C}^{\mathbb{Z}}&\longrightarrow&\mathcal{O}_{d},\\\
(a_{n})&\longmapsto&\mathfrak{B}\displaystyle\sum_{n=-\infty}^{\infty}a_{n}t^{n},\end{array}$

where the $\mathfrak{B}\sum_{n=-\infty}^{\infty}a_{n}t^{n}$ denotes the Borel-
resummation of $\sum_{n=-\infty}^{\infty}a_{n}t^{n}$, is left
$\mathcal{A}$-linear.

###### Proof.

Since the Borel and Laplace transforms are both left $\mathcal{A}$-linear, the
Borel resummation is left $\mathcal{A}$-linear also. ∎

######  Example 2.9.

Let $\mathcal{O}_{2}$ be a space of analytic functions in two variables. Then
$\mathcal{O}_{2}$ becomes a left $\mathcal{A}_{2}$-module by

$\begin{array}[]{ll}\partial_{1}f(x,\,t)=f(x+1,\,t)-f(x,\,t),&X_{1}f(x,\,t)=xf(x-1,\,t);\\\
\partial_{2}f(x,\,t)=f_{t}(x,\,t),&X_{2}f(x,\,t)=tf(x,\,t).\end{array}$ (2.4)

Such a left $\mathcal{A}_{2}$-module is denoted by $\mathcal{O}_{\Delta d}$.

######  Remark 2.10.

In addition to the above notation $\mathcal{O}_{\Delta d}$, we shall liberally
adopt notations such as $\mathcal{O}_{dd}$, $\mathcal{O}_{d\Delta}$,
$\mathcal{O}_{\Delta\Delta}$ to represent different left
$\mathcal{A}_{2}$-module structures endowed on $\mathcal{O}_{2}$ with respect
to the corresponding operators.

We next extend the idea of space of bilateral sequences of numbers
$\mathbb{C}^{\mathbb{Z}}$ in Example 2.5 to the space of bilateral sequences
of analytic functions $\mathcal{O}^{\mathbb{Z}}$.

######  Example 2.11.

Let $\mathcal{O}^{\mathbb{Z}}$ be the space of all bilateral sequences of
analytic functions with appropriately imposed growth restriction. Let

$(f_{n})=(\cdots,\,f_{-1},\,f_{0},\,f_{1},\,\cdots)$

denote a bilateral infinite sequence of analytic functions. Then
$\mathcal{O}^{\mathbb{Z}}$ becomes a left $\mathcal{A}_{2}$-module by

$\begin{array}[]{ll}(\partial_{1}f)_{n}(x)=f^{\prime}_{n}(x),&(X_{1}f)_{n}(x)=xf_{n}(x)\\\
(\partial_{2}f)_{n}(x)=(n+1)f_{n+1}(x),&(X_{2}f)_{n}(x)=f_{n-1}(x),\end{array}$

for all $n$ and for all $x$. This left $\mathcal{A}_{2}$-module is denoted by
$\mathcal{O}_{d}^{\mathbb{Z}}$.

######  Example 2.12.

Let $\mathcal{O}_{\Delta}^{\mathbb{Z}}$ be the space of all bilateral
sequences of analytic functions with an appropriately imposed growth
restriction. Let

$(f_{n})=(\cdots,\,f_{-1},\,f_{0},\,f_{1},\,\cdots)$

denote a bilateral infinite sequence of $f_{k}\in\mathcal{O}$. Then
$\mathcal{O}^{\mathbb{Z}}$ becomes a left $\mathcal{A}_{2}$-module, by

$\begin{array}[]{ll}(\partial_{1}f)_{n}(x)=\Delta
f_{n}(x)=f_{n}(x+1)-f_{n}(x),&(X_{1}f)_{n}(x)=xf_{n}(x-1)\\\
(\partial_{2}f)_{n}(x)=(n+1)f_{n+1}(x),&(X_{2}f)_{n}(x)=f_{n-1}(x),\end{array}$
(2.5)

for all bilateral sequences of analytic functions $(f_{n})$. This left
$\mathcal{A}_{2}$-module is denoted by $\mathcal{O}_{\Delta}^{\mathbb{Z}}$.

Similar to Theorem 2.6, we have the following:

######  Theorem 2.13 ((Generating function of sequences in
$\mathcal{O}_{d}$)).

The following z-transform

$\begin{array}[]{rcl}\mathfrak{z}:\mathcal{O}^{\mathbb{Z}}_{d}&\longrightarrow&\mathcal{O}_{dd},\\\
(f_{n})&\longmapsto&\mathfrak{B}\displaystyle\sum_{n=-\infty}^{\infty}f_{n}(x)\,t^{n}\end{array}$
(2.6)

is left $\mathcal{A}_{2}$-linear.

######  Theorem 2.14 ((Generating function of sequences in
$\mathcal{O}_{\Delta}$)).

The following z-transform

$\begin{array}[]{rcl}\mathfrak{z}_{\triangle}:\mathcal{O}^{\mathbb{Z}}_{\Delta}&\longrightarrow&\mathcal{O}_{\Delta
d},\\\
(f_{n})&\longmapsto&\mathfrak{B}\displaystyle\sum_{n=-\infty}^{\infty}f_{n}(x)\,t^{n}\end{array}$
(2.7)

is left $\mathcal{A}_{2}$-linear.

One needs to make a slight modification of defining rules for
$\mathcal{A}$-module in order to study one-sided (ordinary) sequences. The
motivation of this rule will be explained in Remark 6.3 after we introduce the
Poisson transform.

######  Example 2.15.

Let $\mathbb{C}^{\mathbb{N}}$ be the space of sequences

$a=(a_{n})=(a_{1},\,a_{2},\,a_{3},\,\cdots)$

of complex numbers. It is equipped with the structure of a left
$\mathcal{A}$-module by defining

$\begin{array}[]{l}(Xa)_{n}=na_{n-1}\\\ (\partial a)_{n}=a_{n+1}\quad\mbox{for
all sequences }(a_{n}).\end{array}$

Similarly, let $\mathcal{O}^{\mathbb{N}_{0}}$ be the space of all sequences of
analytic functions with appropriately imposed growth restriction. Let

$(f_{n})=(f_{0},\,f_{1},\,f_{2},\,\cdots)$

denote an infinite sequence of $f_{k}\in\mathcal{O}$. Then
$\mathcal{O}^{\mathbb{N}_{0}}$ becomes a left $\mathcal{A}_{2}$-module by

$\begin{array}[]{ll}(\partial_{1}f)_{n}(x)=f^{\prime}_{n}(x),&(X_{1}f)_{n}(x)=xf_{n}(x)\\\
(\partial_{2}f)_{n}(x)=f_{n+1}(x),&(X_{2}f)_{n}(x)=nf_{n-1}(x),\end{array}$
(2.8)

for all sequences of analytic functions $(f_{n})$. This left
$\mathcal{A}_{2}$-module is usually denoted as
$\mathcal{O}_{d}^{\mathbb{N}_{0}}$.

######  Example 2.16.

Let $\mathcal{O}^{\mathbb{N}_{0}}$ be the space of all sequences of analytic
functions with appropriately imposed growth restriction. Let

$(f_{n})=(f_{0},\,f_{1},\,f_{2},\,\cdots)$

denote an infinite sequence of $f_{k}\in\mathcal{O}$. Then
$\mathcal{O}^{\mathbb{N}_{0}}$ becomes a left $\mathcal{A}_{2}$-module by

$\begin{array}[]{ll}(\partial_{1}f)_{n}(x)=f_{n}(x+1)-f_{n}(x),&(X_{1}f)_{n}(x)=xf_{n}(x-1)\\\
(\partial_{2}f)_{n}(x)=f_{n+1}(x),&(X_{2}f)_{n}(x)=nf_{n-1}(x),\end{array}$
(2.9)

for all sequences of analytic functions $(f_{n})$. This left
$\mathcal{A}_{2}$-module is usually denoted as
$\mathcal{O}_{\Delta}^{\mathbb{N}_{0}}$.

We first demonstrate how we use the theory set up so far to derive a classical
generating function for Bessel functions from the Bessel module
$\mathcal{B}_{\nu}$ defined in Definition 3.38 which is isomorphic to
Definition 3.15.

A table summarising all the $D$-module structures relevant to this article is
included in the Appendix C.1.

###  2.2. Integrability criteria

This section deals with the idea of holonomic modules, which can be found in
[17]. A quick revision is included here so that this article is self-
contained. For each multi-index
$\alpha=(\alpha_{1},\ldots,\alpha_{n})\in\mathbb{N}^{n}$,
$X_{1}^{\alpha_{1}}\cdots X_{n}^{\alpha_{n}}$ will be abbreviated as
$X^{\alpha}$. A similar abbreviation is adopted for the $\partial$’s. It is
also a standard notation that
$|\alpha|=|(\alpha_{1},\ldots,\alpha_{n})|=\sum_{j}\alpha_{j}$.

######  Definition 2.17.

The Bernstein filtration of $\mathcal{A}_{n}$ is the filtration
$F^{0}\mathcal{A}_{n}\subset
F^{1}\mathcal{A}_{1}\subset\cdots\subset\mathcal{A}_{n}$ defined by

$F^{k}\mathcal{A}_{n}=\mbox{the }\mathbb{C}\mbox{-vector space spanned by
}\\{X^{\alpha}\partial^{\beta}:|\alpha|+|\beta|\leq k\\}.$

######  Remark 2.18.

The Bernstein filtration enjoys the following properties.

  * •

$F^{i}\mathcal{A}_{n}\cdot F^{j}\mathcal{A}_{n}\subset
F^{i+j}\mathcal{A}_{n}$,

  * •

$[F^{i}\mathcal{A}_{n},F^{j}\mathcal{A}_{n}]\subset F^{i+j-1}\mathcal{A}_{n}$,

  * •

the $k^{\footnotesize\mbox{th}}-$graded piece
$\mathrm{Gr}^{k}\mathcal{A}_{n}=F^{k}\mathcal{A}_{n}/F^{k+1}\mathcal{A}_{n}$
has finite $\mathbb{C}$-dimension for all $k$.

In fact, in the following treatment, the Bernstein filtration can be replaced
by any other filtration having the these properties.

######  Lemma 2.19.

The graded ring of the Bernstein filtration is a commutative
$\mathbb{C}$-algebra. Indeed,

$\mathrm{Gr}\,\mathcal{A}_{n}=\bigoplus_{k}\mathrm{Gr}^{k}\mathcal{A}_{n}\cong\mathbb{C}[X_{1},\ldots,X_{n},\partial_{1},\ldots,\partial_{n}].$

######  Definition 2.20.

Let $M$ be a left $\mathcal{A}_{n}$-module. A good filtration of $M$ is a
filtration $\Gamma$ which satisfies

  * •

$F^{i}\mathcal{A}_{n}\cdot\Gamma^{j}M\subset\Gamma^{i+j}M$,

  * •

there exists $J\in\mathbb{N}$ such that
$F^{i}\mathcal{A}_{n}\cdot\Gamma^{j}M=\Gamma^{i+j}M$ for all $i$, $j>J$.

######  Lemma 2.21.

Let $M$ be a finitely generated left $\mathcal{A}_{n}$-module and $\Gamma$ be
a good filtration of $M$. Then the graded ring of $\Gamma$, that is,

$\mathrm{Gr}\,M=\bigoplus_{k}{\Gamma^{k}M/\Gamma^{k+1}M},$

is a commutative algebra over the commutative ring
$\mathrm{Gr}\,\mathcal{A}_{n}$.

######  Theorem 2.22 ((Hilbert)).

Let $M$ be a finitely generated left $\mathcal{A}_{n}$-module equipped with a
good filtration $\Gamma$. Then there exists a polynomial $P$ such that

$P(k)=\dim_{\mathbb{C}}\Gamma^{k}M\quad\mbox{for every sufficiently large }k.$

######  Definition 2.23.

The polynomial $P$ in the last theorem is called the Hilbert polynomial of
$M$. (It turns out that $P$ is independent of the choice of the good
filtration chosen.). If $P$ is of degree $d$ and its leading coefficient is
$m/d!$, then we call $d$ the dimension of $M$ and $m$ the multiplicity of $M$.

######  Example 2.24.

Consider $\mathcal{A}_{n}$ as a left $\mathcal{A}_{n}$-module and the
Bernstein filtration is a good one. Then,

$\dim_{\mathbb{C}}\Gamma^{k}\mathcal{A}_{n}=\binom{k+2n}{2n}$

so that the dimension of $\mathcal{A}_{n}$ is $2n$ and its multiplicity is
$1$.

######  Example 2.25.

Consider

$M=\mathcal{A}_{n}/\mathcal{A}_{n}\partial_{1}+\cdots+\mathcal{A}_{n}\partial_{n},$

which inherits a good filtration $\Gamma$ from the Bernstein filtration. Then,

$\dim_{\mathbb{C}}\Gamma^{k}M=\binom{k+n}{n}$

so that the dimension of $M$ is $n$ and its multiplicity is $1$.

######  Example 2.26.

Consider

$M=\mathcal{A}_{2}/\mathcal{A}_{2}\partial_{1}^{2}+\mathcal{A}_{2}\partial_{2},$

which inherits a good filtration $\Gamma$ from the Bernstein filtration. Then,

$\dim_{\mathbb{C}}\Gamma^{k}M=(k+1)^{2}$

so that the dimension of $M$ and its multiplicity are both equal to $2$.

######  Theorem 2.27 ((Bernstein inequality)).

For each non-trivial finitely generated left $\mathcal{A}_{n}$-module $M$, its
dimension $d(M)$ satisfies

$n\leq d(M)\leq 2n.$

######  Definition 2.28.

A finitely generated left $\mathcal{A}_{n}$-module is holonomic if it is
nontrivial and its dimension is $n$.

######  Example 2.29.

$\mathcal{A}_{n}$ is not a holonomic left $\mathcal{A}_{n}$-module. However,
$M=\mathcal{A}_{n}/(\mathcal{A}_{n}\partial_{1}+\cdots+\mathcal{A}_{n}\partial_{n})$
is a holonomic left $\mathcal{A}_{n}$-module.

######  Example 2.30 ((Integrability)).

Let $p(X_{1},X_{2})$, $q(X_{1},X_{2})\in\mathbb{C}[X_{1},X_{2}]$ be
polynomials in two variables. Denote their formal partial derivatives by
subscripts. Consider the left $\mathcal{A}_{2}$-module

$M=\dfrac{\mathcal{A}_{2}}{\mathcal{A}_{2}(\partial_{1}+p(X_{1},X_{2}))+\mathcal{A}_{2}(\partial_{2}+q(X_{1},X_{2}))}.$

  * •

If $p_{2}(X_{1},X_{2})\neq q_{1}(X_{1},X_{2})$, then $M$ is trivial and is
therefore not holonomic.

  * •

If $p_{2}(X_{1},X_{2})=q_{1}(X_{1},X_{2})$, then $M$ can be analyzed in a way
similar to the previous examples, and it turns out that $M$ is holonomic with
its (Hilbert) dimension two and multiplicity one in this case as described in
Example 2.25.

The condition $p_{2}(X_{1},X_{2})=q_{1}(X_{1},X_{2})$ is known as an
integrability condition. Integrability conditions give rise to holonomic
$\mathcal{A}_{n}$-modules, which are important because their “solution space”
(see below) has finite $\mathbb{C}$-dimension. Such a $D$-module thus
ensembles the basic features of an ODE. The examples of $D$-modules which we
study in this article are going to be holonomic.

###  2.3. Banach algebra and Poincaré-Perron theory

In order to study the convergence issue of Bessel’s generating functions that
are results of the $z$-transforms introduced in the last subsection, we extend
the classical Poincaré-Perron theory from sequences of complex numbers to
sequences in Banach algebra [69]. In our applications we only need to consider
sequences in Banach algebra satisfying second order linear difference
equations. However, our adaption can clearly be extended to higher order
difference equations.

######  Theorem 2.31 ((Poincaré-Perron theorem in Banach algebras)).

Let ($\mathbb{B},\|\cdot\|$) be a Banach algebra. Suppose that the $a$,
$b\in\mathbb{B}$ and $(p_{n})$, $(q_{n})$ are sequences in $\mathbb{B}$
converging to zero as $n\to\infty$. If the characteristic equation of the
associated difference equation

$y_{n+2}+ay_{n+1}+by_{n}=0,$ (2.10)

of the difference equation

$y_{n+2}+(a+p_{n})y_{n+1}+(b+q_{n})y_{n}=0,$ (2.11)

has roots $\lambda_{1},\,\lambda_{2}$ distinct in norms, then for each
solution $y_{n}$ to (2.11), the limit

$\lim_{n\to\infty}\frac{\|y_{n+1}\|}{\|y_{n}\|}=\|\lambda_{j}\|,$

holds, for some $j=1,\,2$.

Proofs of Poincaré’s original statement can be found from either Gel’fond [29,
§5.5.2] or Milne-Thomson [48, XVII]. The proof of the above theorem is a
straightforward adaptation in the proof of Poincaré’s original statement with
$\mathbb{B}$ in place of $\mathbb{C}$. Since the proof of Poincaré’s theorem
in full generality is long and the its adaption in Banach algebra $\mathbb{B}$
is almost verbatim except for the aforementioned replacement of “$|\cdot|$” by
“$\|\cdot\|$” in all the computation and estimates, one can verify that a
shorter proof of Poincaré’s theorem available for third order linear
difference equations discussed in [48] can easily be adopted to the Banach
algebra $\mathbb{B}$ setting. Indeed, there are different adaptations of
Poincaré-Perron theory in the literature. For example, the author in [57]
considered objects that assume values in certain abelian groups in place of
complex numbers. Although the asymptotic of the classical Bessel function
$J_{\nu+n}$ as $n\to\infty$ is well-known, a corresponding result for the
difference Bessel function $J^{\Delta}_{\nu+n}$ is unknown and appears highly
non-trivial. The adaption of Poincaré’s theorem in a Banach algebra setting
can easily produce the needed, though rough, asymptotic for the generating
function of difference Bessel functions in to be stated in Theorem 3.39.

##  3\. Bessel modules and Bessel functions

In this section, we investigate the generating functions of the Bessel
functions. We will construct an important holonomic left
$\mathcal{A}_{2}$-module, called the Bessel module, from which we obtain a
pair of PDEs that the generating function satisfy.

######  Definition 3.1 ([13]).

Let $M,N$ be left $\mathcal{A}_{n}$-modules. A solution of $M$ in $N$ is a
left $\mathcal{A}_{n}$-linear maps from $M$ to $N$. In particular, the set of
all solutions of $M$ in $N$ is denoted by

$\mathop{\rm Hom}(M,N).$

######  Theorem 3.2 ([13]).

Let $\mathcal{D}=\langle X\partial,X\rangle$ be the subalgebra of
$\mathcal{A}$ generated by $X\partial$ and $X$. Then there exists an element
$S=\sum_{k=0}^{\infty}a_{k}X^{k}$ in the $X$-adic completion $\mathbb{C}[[X]]$
so that the map

$\mathcal{D}/\mathcal{D}L\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\mathcal{D}/\mathcal{D}(X\partial-\lambda)$

is left $\mathcal{A}$-linear, i.e., it gives a solution of
$\mathcal{D}/\mathcal{D}L$ in $\mathcal{D}/\mathcal{D}(X\partial-\lambda)$.

The element $S$ in the above theorem is a Weyl-algebraic analogue of a
Frobenius series expansion of a solution of $L$. In fact, each left
$\mathcal{A}$-linear map
$\mathcal{A}/\mathcal{A}L\longrightarrow\mathcal{A}/\mathcal{A}K$ induces a
map between solution sets

$\mathop{\rm Hom}(\mathcal{A}/\mathcal{A}K,N)\longrightarrow\mathop{\rm
Hom}(\mathcal{A}/\mathcal{A}L,N).$

Thus finding solutions of a differential operator $L$ in terms of solutions of
another differential operator $K$ is equivalent to identifying an
$\mathcal{A}$-linear map

$\mathcal{A}/\mathcal{A}L\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\mathcal{A}/\mathcal{A}K.$

For instance, finding “power series solutions” of a differential operator $L$
is equivalent to a left $\mathcal{A}$-linear map
$\mathcal{A}/\mathcal{A}L\to\mathcal{A}/\mathcal{A}\partial$; finding a
“Frobenius series solution” of $L$ is equivalent to a left
$\mathcal{A}$-linear map
$\mathcal{A}/\mathcal{A}L\to\mathcal{A}/\mathcal{A}(X\partial-\nu)$; and so
on. We refer to [13] for the details. We illustrate the above discussion with
the following examples that will be needed for the rest of this paper.

###  3.1. Weyl series solutions

######  Example 3.3 ((Weyl-exponential)).

Let $a\in\mathbb{C}$. A solution of $L:=\partial-a$, i.e. a solution of
$\mathcal{A}/\mathcal{A}(\partial-a)$ in any other $D$-module, is called a
Weyl-exponential. In this example, we consider solutions of
$\mathcal{A}/\mathcal{A}(\partial-a)$ in $\mathcal{A}/\mathcal{A}\partial$,
that is, we look for solutions as “power series” of $X$, which are called Weyl
power series to the Weyl exponential. For this purpose we compute the element
$S=\sum_{k}c_{k}X^{k}$ which gives rise to a left $\mathcal{A}$-linear map

$\mathcal{A}/\mathcal{A}(\partial-a)\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\overline{\mathcal{A}/\mathcal{A}\partial}$

guaranteed by the study in [13, §4.1]. For each $k\geq 1$, since $\partial
X^{k}=X^{k}\partial+kX^{k-1}$, we have that

$\begin{split}0=(\partial-a)\left(\sum_{k=0}^{\infty}c_{k}X^{k}\right)&=\sum_{k=0}^{\infty}c_{k}X^{k}\partial+\sum_{k=1}^{\infty}c_{k}kX^{k-1}-a\sum_{k=0}^{\infty}c_{k}X^{k}\\\
&=\sum_{k=0}^{\infty}(c_{k+1}(k+1)-ac_{k})X^{k}.\quad\qquad\mod\mathcal{A}\partial\end{split}$

The above formula gives recurrence of $c_{k}$ from the relation

$kc_{k}-ac_{k-1}=0,\qquad k\geq 1.$

By choosing $c_{0}=1,$ we obtain $c_{k}=a^{k}/k!$. Thus we denote

$S(X)={\mathop{\rm E}}(aX):=\sum_{k=0}^{\infty}\frac{a^{k}}{k!}X^{k}$ (3.1)

to be the Weyl exponential series.

######  Example 3.4.

(Weyl Sine and Cosine) The Weyl-sine and the Weyl-cosine are solutions of
$L:=\partial^{2}+1$, i.e., solutions of the $D$-module
$\mathcal{A}/\mathcal{A}(\partial^{2}+1)$. To express them as “Weyl power
series”, i.e., in terms of powers of $X$, we again let $K=\partial$, and
compute the element $S=\sum c_{k}X^{k}$ which yields the left $\mathcal{A}$

$\mathcal{A}/\mathcal{A}(\partial^{2}+1)\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\overline{\mathcal{A}/\mathcal{A}\partial},$ (3.2)

which are solutions of $\mathcal{A}/\mathcal{A}(\partial^{2}+1)$ in
$\mathcal{A}/\mathcal{A}\partial$. Since
$\partial^{2}+1=(\partial-i)(\partial+i)$, we apply Example 3.3 with the
choice $a=-i$ from the Example 3.3. Then it follows immediately that

$(\partial-i)(\partial+i){\mathop{\rm E}}(-iX)=0\qquad\mod\partial$

and since $\partial^{2}+1=(\partial+i)(\partial-i)$ so that

$(\partial+i)(\partial-i){\mathop{\rm E}}(iX)=0\qquad\mod\partial,$

where we have chosen $a=i$ from the Example 3.3. Thus both ${\mathop{\rm
E}}(iX)$ and ${\mathop{\rm E}}(-iX)$ in $\mathcal{A}$ are possible choices of
$S$ in (3.2). Hence both

$\mathop{\rm Sin}(X):=\frac{1}{2i}\big{(}{\mathop{\rm E}}(iX)-{\mathop{\rm
E}}(-iX)\big{)}$ (3.3)

and

$\mathop{\rm Cos}(X):=\frac{1}{2}\big{(}{\mathop{\rm E}}(iX)+{\mathop{\rm
E}}(-iX)\big{)}$ (3.4)

are different choices of elements of $\mathcal{A}$, called the Weyl-Sine and
Cosine series, that can also serve for the $S$ in (3.2). corresponding

Alternatively, we may compute the Weyl Sine and Cosine series $S=\sum
c_{k}X^{k}$ directly by

$\begin{split}0&=(\partial-i)(\partial+i)\sum_{k=0}^{\infty}c_{k}X^{k}\\\
&=\sum_{k=0}^{\infty}c_{k}X^{k}\partial^{2}+2\sum_{k=1}^{\infty}c_{k}kX^{k-1}\partial+\sum_{k=2}^{\infty}c_{k}k(k-1)X^{k-2}+\sum_{k=0}^{\infty}c_{k}X^{k}\\\
&=\sum_{k=0}^{\infty}c_{k+2}(k+1)(k+2)X^{k}+\sum_{k=0}^{\infty}c_{k}X^{k}.\qquad\mod\partial\end{split}$

That is,

$c_{k+2}=\frac{-c_{k}}{(k+1)(k+2)},\quad k\geq 0.$

Thus, if $c_{0}=0,c_{1}=1$, we obtain Weyl Sine series

$\mathop{\rm
Sin}X=X-\frac{X^{3}}{3!}+\frac{X^{5}}{5!}-\cdots+\frac{(-1)^{k+1}}{(2k-1)!}X^{2k-1}+\cdots.$

If $c_{0}=1,c_{1}=0$, we obtain Weyl (power) Cosine series

$\mathop{\rm
Cos}X=1-\frac{X^{2}}{2!}+\frac{X^{4}}{4!}-\cdots+\frac{(-1)^{k}}{(2k)!}X^{2k}+\cdots.$

Hence they serve as $S$ with “appropriately imposed initial conditions”.

We now introduce one of the main subjects of study in this paper.

######  Example 3.5 ((Weyl Bessel of order $\nu$)).

Let $\nu\in\mathbb{C}$ such that $2\nu$ is not an integer, and let

$L=(X\partial)^{2}+X^{2}-\nu^{2}.$ (3.5)

Let $\mathcal{D}=\langle X\partial,X\rangle$ be the subalgebra of
$\mathcal{A}$ generated by $X\partial$ and $X$. Then the Weyl Bessel of order
$\nu$ is a solution of $L$ in a left $\mathcal{D}$-module $N$ , i.e., a
solution of $\mathcal{D}/\mathcal{D}L$ in $N$. To find the “Frobenius series
expansion” of this Weyl Bessel, we look for a solution of
$\mathcal{D}/\mathcal{D}L$ in the $D$-module
$\overline{\mathcal{D}/\mathcal{D}(X\partial-\nu)}$, i.e., we find the element
$S=\sum_{k=0}^{\infty}c_{k}X^{k}$ that gives the left $\mathcal{A}$-linear map

$\mathcal{D}/\mathcal{D}L\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\overline{\mathcal{D}/\mathcal{D}(X\partial-\nu)}.$
(3.6)

Since

$(X\partial+\nu)(X\partial-\nu)X^{k}=X^{k}[(X\partial+\nu)+2k](X\partial-\nu)+(2\nu
k+k^{2})X^{k},$

holds for every $k\geq 0$ (as well as trivially when $k=0$), we have

$\begin{split}0&=L\sum_{k=0}^{\infty}c_{k}X^{k}=\sum_{k=0}^{\infty}c_{k}[(X\partial+\nu)(X\partial-\nu)+X^{2}]X^{k}\\\
&=\sum_{k=0}^{\infty}c_{k}X^{k}[(X\partial+\nu)+2k](X\partial-\nu)+\sum_{k=1}^{\infty}{c_{k}k(2\nu+k)X^{k}}+\sum_{k=0}^{\infty}c_{k}X^{k+2}\\\
&=\sum_{k=1}^{\infty}c_{k}k(2\nu+k)X^{k}+\sum_{k=0}^{\infty}c_{k}X^{k+2}\hfill\mod
X\partial-\nu\\\
&=c_{1}(1+2\nu)X+\sum_{k=2}^{\infty}(c_{k}k(2\nu+k)+c_{k-2})X^{k},\hfill\mod
X\partial-\nu\end{split}$

which implies that $c_{1}=c_{3}=c_{5}=\cdots=0$, and by choosing
$c_{0}\Gamma(\nu+1)=2^{-\nu},$ we must have

$c_{2k}=\frac{(-1)^{k}}{2^{\nu+2k}k!\Gamma(\nu+k+1)},\quad k\geq 0.$

This shows that

$S=\sum_{k=0}^{\infty}\frac{(-1)^{k}X^{2k}}{2^{\nu+2k}k!\Gamma(\nu+k+1)}$

can be chosen as such a Weyl series $S$ in (3.6) that serves as a solution for
(3.6). We denote this specific Weyl Bessel series by $\mathcal{J}_{\nu}$ as
solution in $\mathcal{D}/\mathcal{D}(X\partial-\nu)$, or more generally, in a
left $\mathcal{A}$-module $N$, which can be considered as the Weyl algebraic
version of the classical Bessel $J_{\nu}$. We caution that the reader not to
confuse the series form of $\mathcal{J}_{\nu}$ in contrast to the original
Bessel function of order $J_{\nu}$ since there is no factor “$X^{\nu}$” that
would appear there unless when the $\nu$ reduces to an integer.

###  3.2. Transmutation formulae

After defining the Weyl Bessel as a solution of
$L=L_{\nu}:=(X\partial)^{2}+X^{2}-\nu^{2}$ in the last subsection, we observe
the following transmutation formulae for $L$. They can be verified directly,
but they do not appear to be found immediately amongst the formulae in
Dereziński and Majewski [20, 21, 22]141414Please see §9 Conclusion for
discussion related to the work of Infeld and Hull [35].. As we shall see that
they are instrumental in the development of Bessel module to be defined in the
next subsection.

######  Proposition 3.6.

For each $\nu\in\mathbb{C}$, we have

$\displaystyle[(X\partial)^{2}+X^{2}-(\nu-1)^{2}]\left(\partial+\dfrac{\nu}{X}\right)$
$\displaystyle=\left(\partial+\dfrac{\nu-2}{X}\right)[(X\partial)^{2}+X^{2}-\nu^{2}],$
(3.7)
$\displaystyle[(X\partial)^{2}+X^{2}-(\nu+1)^{2}]\left(\partial-\dfrac{\nu}{X}\right)$
$\displaystyle=\left(\partial-\dfrac{\nu+2}{X}\right)[(X\partial)^{2}+X^{2}-\nu^{2}],$
(3.8)

which induce the following left $\mathcal{A}$-linear maps

$\begin{array}[]{rcl}\mathcal{A}/\mathcal{A}((X\partial)^{2}+X^{2}-(\nu-1)^{2})&\stackrel{{\scriptstyle\times(\partial+\frac{\nu}{X})}}{{\longrightarrow}}&\mathcal{A}/\mathcal{A}((X\partial)^{2}+X^{2}-\nu^{2}),\\\
\mathcal{A}/\mathcal{A}((X\partial)^{2}+X^{2}-(\nu+1)^{2})&\stackrel{{\scriptstyle\times(\partial-\frac{\nu}{X})}}{{\longrightarrow}}&\mathcal{A}/\mathcal{A}((X\partial)^{2}+X^{2}-\nu^{2}).\end{array}$

######  Corollary 3.7.

Let $\nu\in\mathbb{C}$, the classical Bessel functions satisfy the following
recurrence relations.

$\displaystyle J^{\prime}_{\nu}(x)+\dfrac{\nu J_{\nu}(x)}{x}$
$\displaystyle=J_{\nu-1}(x),$ (3.9) $\displaystyle
J^{\prime}_{\nu}(x)-\dfrac{\nu J_{\nu}(x)}{x}$ $\displaystyle=-J_{\nu+1}(x).$
(3.10)

###### Proof.

This proof serves to provide an alternative to the classical proof found in
most literature. It is a direct consequence from the previous Proposition with
the $D$-module structure on $\mathcal{O}$ introduced from the Example 2.3. Let
$J_{\nu}$ be the classical Bessel function to $L_{\nu}$ with the asymptotic
behaviour $\lim_{x\to 0^{+}}{J_{\nu}(x)}/{x^{\nu}}={1}/{2^{\nu}\Gamma(\nu+1)}$
[64, p. 43]. Then one applies the both sides of (3.7). It implies that
$(\partial+\nu/X)J_{\nu}=CJ_{\nu-1}+DJ_{-\nu+1}$ holds for some constants
$C,D$. One can easily verify that the $C=1,\ D=0$ from the asymptotic
$J^{\prime}_{\nu}(x)/x^{\nu-1}\sim 1/2^{\nu-1}\Gamma(\nu)$ [64, p. 43]. This
verifies the (3.9). The verification of (3.10) is similar.

∎

######  Corollary 3.8.

Let $\nu\in\mathbb{C}$ and a left $\mathcal{D}$-module $N$ be given. Then
there exists a sequence $\mathrm{(}\mathfrak{J}_{\nu+n}\mathrm{)}$ in
$N^{\mathbb{Z}}$ and for each $n$ the $\mathfrak{J}_{\nu+n}$ is a solution to
$L_{n}=(X\partial)^{2}+X^{2}-(\nu+n)^{2}$ in $N$, such that we have

$\displaystyle
X\partial\,\mathfrak{J}_{\nu+n}+(\nu+n)\mathfrak{J}_{\nu+n}-X\mathfrak{J}_{\nu+n-1}=0,$
(3.11) $\displaystyle
X\partial\,\mathfrak{J}_{\nu+n}-(\nu+n)\mathfrak{J}_{\nu+n}+X\mathfrak{J}_{\nu+n+1}=0.$
(3.12)

######  Remark 3.9.

Reader are cautioned not to confuse the notation $\mathfrak{J}_{\nu}$
introduced in the last Corollary with notation $\mathcal{J}_{\nu}$ that
denotes the Weyl-Bessel series introduced in Example 3.5 that is defined as a
solution with respect to the (3.6).

###### Proof.

Recall that $L_{\nu}=(X\partial)^{2}+X^{2}-\nu^{2}$. Let us replace $\nu$ by
$\nu+n$ in (3.7). Then it is clear that the right-side of (3.7) annihilates
$\mathfrak{J}_{\nu+n}$. However the expression
$(\partial+\nu/X)\mathfrak{J}_{\nu+n}=(X\partial+\nu)/X\,(X\mathfrak{J}_{\nu}$)
on the left-side must be annihilated by $L_{\nu+n-1}$. Hence
$X\partial\,\mathfrak{J}_{\nu+n}+(\nu+n)\mathfrak{J}_{\nu+n}=\pm
X\mathfrak{J}_{\nu+n-1}$. This gives the (3.11) except for the “$\mp$” in
front of the third term $X\mathfrak{J}_{\nu+n-1}$. We have chosen the minus
sign over the plus sign in accordance with the corresponding formula for the
classical Bessel function $J_{\nu}$. This establishes the (3.11). The proof
for (3.12) is similar. ∎

######  Corollary 3.10 ((Three-term recurrence)).

Let $\nu\in\mathbb{C}$ and a left $\mathcal{D}$-module $N$ be given. Then
there exists a sequence $\mathrm{(}\mathfrak{J}_{\nu+n}\mathrm{)}$ such that
for each $n$ the $\mathfrak{J}_{\nu+n}$ is a solution to
$L_{n}=(X\partial)^{2}+X^{2}-(\nu+n)^{2}$ in $N$, such that we have

$2(\nu+n)\,\mathfrak{J}_{\nu+n}-X\mathfrak{J}_{\nu+n-1}-X\mathfrak{J}_{\nu+n+1}=0$
(3.13)

as an element in $N$.

###  3.3. Bessel modules

Recall that the Weyl algebra $\mathcal{A}_{2}$ is the free
$\mathbb{C}$-algebra
$\mathbb{C}\langle\partial_{1},\partial_{2},X_{1},X_{2}\rangle$ subject to

$\partial_{i}\partial_{j}-\partial_{j}\partial_{i}=0,\quad
X_{i}X_{j}-X_{j}X_{i}=0,\quad\partial_{i}X_{j}-X_{j}\partial_{i}=\delta_{ij},$

for $i,j=1,2$. Then we define

######  Definition 3.11.

We denote the algebra obtained by adjoining the new elements
$1/X_{1},\,1/X_{2}$ to the Weyl algebra $\mathcal{A}_{2}$ by

$\mathcal{A}_{2}(1/X_{1},\,1/X_{2}):=\mathbb{C}\langle\partial_{1},\partial_{2},X_{1},\dfrac{1}{X_{1}},X_{2},\dfrac{1}{X_{2}}\rangle.$

Let $\nu\in\mathbb{C}$. Motivated by Corollary 3.8, we now consider the left
$\mathcal{A}_{2}$-module generanted by the following elements

$X_{1}\partial_{1}+(X_{2}\partial_{2}+\nu)-X_{1}X_{2},\qquad
X_{1}\partial_{1}-(X_{2}\partial_{2}+\nu)+\frac{X_{1}}{X_{2}}.$ (3.14)

With the (3.14) in mind, we define

######  Definition 3.12 ((Bessel module of order $\nu$)).

Let $\nu\in\mathbb{C}$. The left $\mathcal{A}_{2}$-module151515The
$\mathcal{B}_{\nu}$ is in fact a quotient module. We have left out the word
“quotient” when no confusion can arise.

$\mathcal{B}_{\nu}=\frac{\mathcal{A}_{2}}{\mathcal{A}_{2}(X_{1}\partial_{1}+(\nu+X_{2}\partial_{2})-X_{1}X_{2})+\mathcal{A}_{2}(X_{1}\partial_{1}-(\nu+X_{2}\partial_{2})+{X_{1}}/{X_{2}})}$
(3.15)

is called the Bessel module of order $\nu$.

First we check that the $\mathcal{B}_{\nu}$ is indeed holonomic.

######  Proposition 3.13.

Let $\nu\in\mathbb{C}$. The Bessel module $\mathcal{B}_{\nu}$ is a left
$\mathcal{A}_{2}$-module 161616We omit the description that
$\mathcal{B}_{\nu}$ is a left $\mathcal{A}_{2}$-quotient module for the sake
of simplicity. with dimension $2$ and multiplicity $1$. In particular, it is
holonomic.

###### Proof.

The two generators from Definition 3.12 that generate the Bessel module
$\mathcal{B}_{\nu}$ satisfy the integrability criterion stated in Example
2.25. Hence the Bessel module has dimension $2$ and multiplicity $1$. ∎

The following theorem shows that one can recover from $\mathcal{B}_{\nu}$ the
ODE satisfied by Bessel functions.

######  Proposition 3.14 ((Bessel ODE module of order $\nu$)).

The map

$\dfrac{\mathcal{A}_{2}}{\mathcal{A}_{2}[(X_{1}\partial_{1})^{2}+X_{1}^{2}-(\nu+X_{2}\partial_{2})^{2}]}\longrightarrow\mathcal{B}_{\nu}$

is a well-defined left $\mathcal{A}_{2}$-linear surjection.

###### Proof.

Recall that $\mathcal{B}_{\nu}$ is generated by the two elements

$X_{1}\partial_{1}+X_{2}\partial_{2}-X_{1}X_{2}+\nu,\qquad
X_{1}\partial_{1}-X_{2}\partial_{2}+\frac{X_{1}}{X_{2}}-\nu$

as in (3.14), so as an element in $\mathcal{B}_{\nu}$, we have

$\displaystyle(X_{1}\partial_{1})^{2}+X_{1}^{2}-(\nu+X_{2}\partial_{2})^{2}$
$\displaystyle=(X_{1}\partial_{1})^{2}+X_{1}^{2}+(\nu+X_{2}\partial_{2})(X_{1}\partial_{1}-X_{1}X_{2})$
$\displaystyle=(X_{1}\partial_{1})(X_{1}\partial_{1}+\nu+X_{2}\partial_{2})+(X_{1}X_{2})(\frac{X_{1}}{X_{2}}-\nu-\partial_{2}X_{2})$
$\displaystyle=(X_{1}\partial_{1})(X_{1}X_{2})+(X_{1}X_{2})(X_{2}\partial_{2}-X_{1}\partial_{1}-\partial_{2}X_{2})$
$\displaystyle=(X_{1}\partial_{1})(X_{1}X_{2})+(X_{1}X_{2})(-\partial_{1}X_{1})$
$\displaystyle=0.$

Hence the map is well defined. ∎

We may rephrase the above result in the following form.

######  Proposition 3.15.

Let $\nu\in\mathbb{C}$ and $M$ be a left $\mathcal{A}$-module. If the sequence
$(f_{n})$ is a solution of $\mathcal{B}_{\nu}$ in $M^{\mathbb{Z}}$, then for
each $n\in\mathbb{Z}$, $f_{n}$ is a solution of
$(X\partial)^{2}+X^{2}-(\nu+n)^{2}$ in $M$. In other words, the element
$(X_{1}\partial_{1})^{2}+X_{1}^{2}-(\nu+X_{2}\partial_{2})^{2}$ is equivalent
to $0$ in $\mathcal{B}_{\nu}$.

###  3.4. Characteristics of Bessel modules

We will see later that the Bessel module $\mathcal{B}_{\nu}$ gives rise to
PDEs satisfied by the generating functions of the Bessel functions
$(J_{\nu+n})$ in §3.5 as well as of the recently discovered difference Bessel
functions $(J^{\Delta}_{\nu+n})$ in §3.6 below. The following theorem is an
algebraic analogue of the method of characteristics, which will help in
solving such PDEs.

######  Theorem 3.16.

Let $\nu\in\mathbb{C}$ and $\mathcal{B}_{\nu}$ be defined above. Then there
exists a left $\mathcal{A}_{2}(1/X_{j=1,2})$-linear map $S$

$\mathcal{B}_{\nu}\xrightarrow{\times
S}\overline{\mathcal{A}_{2}/\big{[}\mathcal{A}_{2}\partial_{1}+\mathcal{A}_{2}(X_{2}\partial_{2}+\nu)\big{]}}.$
(3.16)

In fact, we can take

$S={\mathop{\rm
E}}\Big{[}\frac{X_{1}}{2}\Big{(}X_{2}-\frac{1}{X_{2}}\Big{)}\Big{]},$ (3.17)

where ${\mathop{\rm E}}$ denotes the Weyl exponential series introduced in
Example 3.1.

###### Proof.

Define new symbols $\Theta_{1},\Theta_{2},Y_{1},Y_{2}$ by

$Y_{1}=X_{1}X_{2},\qquad Y_{2}=\frac{X_{1}}{X_{2}},\qquad
2\Theta_{1}=\frac{\partial_{1}}{X_{2}}+\frac{\partial_{2}}{X_{1}},\qquad
2\Theta_{2}=X_{2}\partial_{1}-\frac{X^{2}_{2}\partial_{2}}{X_{1}}.$

Direct verification yields

$[\Theta_{1},Y_{1}]=1,\quad[\Theta_{2},Y_{2}]=1,\quad[\Theta_{1},\Theta_{2}]=0,\quad[Y_{1},Y_{2}]=0,\quad[\Theta_{i},\,Y_{j}]=0,\
i\not=j.$ (3.18)

The elements in (3.14) become

$2Y_{1}\Theta_{1}-Y_{1}+{\nu}\qquad\mbox{and}\qquad
2Y_{2}\Theta_{2}+Y_{2}-{\nu}.$

It follows from the Weyl exponential introduced in Example 3.3 for $a=\nu/2$
that

$\begin{split}(2Y_{1}\Theta_{1}-Y_{1}+\nu)\,{\mathop{\rm
E}}(Y_{1}/2)&=2Y_{1}{\mathop{\rm E}}(Y_{1}/2)\Theta_{1}+\nu{\mathop{\rm
E}}(Y_{1}/2)\\\ &={\mathop{\rm E}}(Y_{1}/2)(2Y_{1}\Theta_{1}+\nu)\\\
&=0.\qquad\qquad\quad\mod\mathcal{A}_{2}(Y_{1}\Theta_{1}+\nu/2)\end{split}$

In general we have $S=F(Y_{2})\cdot{\mathop{\rm E}}(Y_{1}/2)$ for some factor
$F(Y_{2})$ depending on $Y_{2}$ only. Next we require $F(Y_{2})$ to satisfy
that

$(2Y_{2}\Theta_{2}+Y_{2}-\nu)\,F(Y_{2})\cdot{\mathop{\rm
E}}(Y_{1}/2)=0.\qquad\mod\mathcal{A}_{2}(Y_{2}\Theta_{2}-\nu/2)$

That is,

$(2Y_{2}\Theta_{2}+Y_{2}-\nu)\,F(Y_{2})=0,\qquad\qquad\mod\mathcal{A}_{2}(Y_{2}\Theta_{2}-\nu/2).$

It follows from a similar consideration as above that $F(Y_{2})=c{\mathop{\rm
E}}(-Y_{2}/2)$ for some complex number $c$. Without loss of generality we may
choose $c=1$ so that the above map becomes

$S:={\mathop{\rm E}}(Y_{1}/2)\cdot{\mathop{\rm
E}}(-Y_{2}/2)=\sum_{k=0}^{\infty}\frac{X_{1}^{k}X_{2}^{k}}{2^{k}k!}\cdot\sum_{\ell=0}^{\infty}\frac{(-1)^{\ell}X_{1}^{\ell}}{X_{2}^{\ell}2^{\ell}\ell!}={\mathop{\rm
E}}\Big{[}\frac{X_{1}}{2}\Big{(}X_{2}-\frac{1}{X_{2}}\Big{)}\Big{]}$ (3.19)

which is the map asserted in (3.16). ∎

######  Remark 3.17.

We may rewrite the product of two series in $S$ in (3.19) to get

$\sum_{k=0}^{\infty}\frac{X^{k}_{1}X^{k}_{2}}{2^{k}k!}\sum_{\ell=0}^{\infty}\frac{(-1)^{l}X^{\ell}_{1}}{X^{\ell}_{2}2^{\ell}!}=\sum_{n=-\infty}^{\infty}\left(\sum_{\ell=0}^{\infty}\frac{(-1)^{\ell}X^{n+2\ell}_{1}}{2^{n+2\ell}(n+\ell)!\ell!}\right)X_{2}^{n}.$

The “inner sum”

$\mathcal{J}_{n}(X)=\sum_{\ell=0}^{\infty}\frac{(-1)^{\ell}}{(n+\ell)!l!}\left(\frac{X}{2}\right)^{n+2\ell}$

agrees with Example 3.5 after composing with the multiplication by $X^{n}$
from $\mathcal{A}/\mathcal{A}(X\partial-n)$ to
$\mathcal{A}/\mathcal{A}\partial$. It gives a “power series” of the Weyl
Bessel of order $n$. Following the arguments used in Watson [64, pp. 15-16]
that

$\mathcal{J}_{-n}(X)=(-1)^{n}\mathcal{J}_{n}(X)$ (3.20)

holds for each integer $n\geq 1$.

###  3.5. Generating function of classical Bessel functions

As was explained by Watson [64, Chap. 2., §2,.1] that the convergence of the
infinite expansion of the generating function

$e^{\frac{x}{2}\big{(}t-\frac{1}{t}\big{)}}=\sum_{n=-\infty}^{\infty}J_{n}(x)\,t^{n}$
(3.21)

of Bessel function is absolute in powers of $t$, and valid for all $x$ and
$t\not=0$. We shall take this expansion as the basis of our study of various
extensions of this generating function for various types of Bessel functions
in this section.

######  Example 3.18 ((Frobenius series of classical Bessel functions of order
$\nu$)).

Let $\nu$ be a complex number such that $2\nu$ is not an integer, and consider
the the Bessel operator $L=(X\partial)^{2}+X^{2}-\nu^{2}$ as in Example 3.5.
The classical Bessel function of order $\nu$, $J_{\nu}(x)$, is a solution of
$\mathcal{D}/\mathcal{D}L$ in the $D$-module $\mathcal{O}_{d}$. To obtain a
Frobenius series expansion of $J_{\nu}(x)$, we recall from Example 3.5 the
element

$S=\mathcal{J}_{\nu}(X)=\sum_{k=0}^{\infty}\frac{(-1)^{k}X^{2k}}{2^{\nu+2k}k!\Gamma(\nu+k+1)}$
(3.22)

gives the composition map of left $\mathcal{A}$-linear maps

$\mathcal{D}/\mathcal{D}L\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\overline{\mathcal{D}/\mathcal{D}(X\partial-\nu)}\stackrel{{\scriptstyle\times
x^{\nu}}}{{\longrightarrow}}\mathcal{O}_{d},$ (3.23)

and so

$J_{\nu}(x)=x^{\nu}\sum_{k=0}^{\infty}\frac{(-1)^{k}x^{2k}}{2^{\nu+2k}k!\Gamma(\nu+k+1)},$
(3.24)

gives raise to the classical Bessel function of order $\nu$.

On the other hand, if we replace the $\nu$ in (3.12) by $-\nu$, and proceed in
a similar manner for an analogous set up, then we would have obtained a
“linearly-independent” solution of $\mathcal{D}/\mathcal{D}L$ in
$\mathcal{O}_{d}$ given by

$\mathcal{D}/\mathcal{D}L\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\overline{\mathcal{D}/\mathcal{D}(X\partial+\nu)}\stackrel{{\scriptstyle\times
x^{-\nu}}}{{\longrightarrow}}\mathcal{O}_{d},$ (3.25)

where

$S=\mathcal{J}_{-\nu}(X)=\sum_{k=0}^{\infty}\frac{(-1)^{k}X^{2k}}{2^{-\nu+2k}k!\Gamma(-\nu+k+1)}.$
(3.26)

The (3.25) gives rise to the Frobenius series expansion

$J_{-\nu}(x)=x^{-\nu}\sum_{k=0}^{\infty}\frac{(-1)^{k}x^{2k}}{2^{-\nu+2k}k!\Gamma(-\nu+k+1)}.$
(3.27)

The following example illustrates how sequences of analytic functions are
realized from $D$-modules.

######  Example 3.19.

Let $\nu\in\mathbb{C}$ and consider the bilateral sequence of classical Bessel
functions $(J_{\nu+n})_{n}$. Then the maps

$\begin{array}[]{rcl}\mathfrak{j}_{+}:\mathcal{B}_{\nu}&\stackrel{{\scriptstyle\times(J_{\nu+n})}}{{\longrightarrow}}&\mathcal{O}^{\mathbb{Z}}\end{array}$
(3.28)

and

$\begin{array}[]{rcl}\mathfrak{j}_{-}:\mathcal{B}_{\nu}&\stackrel{{\scriptstyle\times(J_{-\nu-n})}}{{\longrightarrow}}&\mathcal{O}^{\mathbb{Z}}\end{array}$
(3.29)

are both left $\mathcal{A}_{2}$-linear.

More generally, a general solution of $\mathcal{B}_{\nu}$ in
$\mathcal{O}^{\mathbb{Z}}$ will be denoted by
$\times\mathrm{(}\mathscr{C}_{\nu+n}\mathrm{)}$. That is,

$\begin{array}[]{rcl}\mathfrak{j}:\mathcal{B}_{\nu}&\stackrel{{\scriptstyle\times(\mathscr{C}_{\nu+n})}}{{\longrightarrow}}&\mathcal{O}^{\mathbb{Z}}\end{array}$
(3.30)

The following proposition and corollary are direct consequences of Corollary
3.8 with $M=\mathcal{O}_{d}$ and $\mathscr{C}_{\nu+n}$.

######  Proposition 3.20 ((Differential-difference formulae)).

Let $\mathcal{B}_{\nu}\ \begin{subarray}{c}\times(\mathcal{C}_{\nu+n})\\\
\longrightarrow\end{subarray}\ \mathcal{O}^{\mathbb{Z}}_{d}$
$\mathrm{(}\mathscr{C}_{\nu+n}\mathrm{)}$ be a solution of $\mathcal{B_{\nu}}$
in $\mathcal{O}^{\mathbb{Z}}_{d}$. Then the
$\mathrm{(}\mathscr{C}_{\nu+n}\mathrm{)}$ satisfies the following classical
differential-difference formulae

$x\mathscr{C}_{\nu}^{\prime}(x)+\nu\mathscr{C}_{\nu}(x)-x\mathscr{C}_{\nu-1}(x)=0,$
(3.31)

and

$x\mathscr{C}^{\prime}_{\nu}(x)-\nu\mathscr{C}_{\nu+n}(x)+x\mathscr{C}_{\nu+1}(x)=0.$
(3.32)

The usage of the notation $\mathscr{C}_{\nu+n}$ follows Sonine’s notation (see
[64, p. 83]) for functions that satisfy the system (3.31) and (3.32).

Here is an application of the proposition to recover two well-known formulae
about the classical Bessel functions. We shall see how this will lead to new
results for difference Bessel functions and half-Bessel modules for difference
reverse Bessel polynomials in §6.2.

One recovers the well-known formulae (3.9) and (3.10) by choosing
$\mathscr{C}_{\nu}=J_{\nu}$ in the above Corollary.

The following classical three-term recursion formula is a consequence of the
above differential-difference formulae (3.31) and (3.32).

######  Corollary 3.21 ((Three-term recurrence)).

Let $\mathcal{B}_{\nu}\ \begin{subarray}{c}\times(\mathcal{C}_{\nu+n})\\\
\longrightarrow\end{subarray}\ \mathcal{O}^{\mathbb{Z}}_{d}$
$\mathrm{(}\mathscr{C}_{\nu+n}\mathrm{)}$ be a solution of $\mathcal{B_{\nu}}$
in $\mathcal{O}^{\mathbb{Z}}_{d}$. Then

$2\nu\mathscr{C}_{\nu}(x)-x\mathscr{C}_{\nu-1}(x)-x\mathscr{C}_{\nu+1}(x)=0.$
(3.33)

Now we come to the main result of this subsection, which is about the
generating function of the classical Bessel functions $(J_{\nu+n}(x))$. It is
a two-variable analytic solution of the holonomic system of PDEs obtained from
the Bessel module $\mathcal{B}_{\nu}$ (3.15).

######  Theorem 3.22 ((Asymptotic behaviour)).

The sequence $(J_{n+\nu})$ is a uniformly 1-Gevrey compacta solution of
$\mathcal{B}_{\nu}$ in $\mathcal{O}_{d}^{\mathbb{Z}}$.

###### Proof.

It follows from the the Corollary 3.21 that $J_{\nu+n}$ satisfies

$y_{n+2}-\dfrac{2(n+\nu+1)}{x}y_{n+1}+y_{n}=0.$

As a result, the modified sequence ($Y_{n}$) where $Y_{n}:=J_{\nu+n}/n!$
satisfies the equation

$Y_{n+2}-\dfrac{2(n+\nu+1)}{(n+2)x}Y_{n+1}+\dfrac{1}{(n+2)(n+1)}Y_{n}=0.$
(3.34)

Now one fixes a compact set $K\subset\mathbb{C}\backslash\\{0\\}$. The space
of functions analytic on $K$ becomes a Banach algebra $\mathbb{B}_{K}$
equipped with the usual sup norm $\|\cdot\|_{K}$.Therefore, we may apply the
generalized Poincaré-Perron theory stated in Theorem 2.31 on the Banach
algebra $\mathbb{B}_{K}$ for the asymptotic behaviour of ($J_{\nu+n}$) on the
compact set $K$. Since the characteristic equation of (3.34) is given by

$\lambda^{2}-2\lambda/x=0,$

so that its roots are $0$ and $2/x$ whose norms are obviously different.
Consequently

$\|J_{n+\nu}/n!\|_{K}=O((2/d(K,0))^{n}).$ (3.35)

That is, according to the Remark 2.8, the sequence $(J_{n+\nu})$ is
$1$-Gevrey. ∎

######  Remark 3.23.

  1. (1)

It is instructive to compare the estimate (3.35) with the well-known
asymptotic

$J_{\nu+n}(x)\sim\frac{e^{-\nu}}{\sqrt{2\pi
n}}\Big{(}\frac{ex}{2n}\Big{)}^{\nu+n}$

as given, see for example, in [28, §5]. An application of Stirling formula to
(3.35) shows that we have the same order of growths from the two asymptotic
estimates. 171717In fact, this can be seen from the elementary growth estimate
[64, p. 225]
$J_{\nu}(x)\approx\exp\big{\\{}\nu+\nu\log(x/2)-(\nu+\frac{1}{2})\log\nu\big{\\}}\cdot\big{[}1/\sqrt{2\pi}+\frac{c_{1}}{\nu}+\frac{c_{2}}{\nu^{2}}+\cdots\big{]}$
holds for a fixed $x$ as $\nu\to\infty$, so that although the “analystic
part”, i.e., the “$\sum_{0}^{\infty}$” part of (3.36) is convergent, while the
“principle part”, or the “$\sum_{-1}^{-\infty}$” is divergent.

  2. (2)

Similarly, one sees that the sequence $(J_{-n+\nu})$ is “uniformly 1-Gevrey
compacta” also. The conclusion remains valid if the standard Bessel functions
are replaced by an arbitrary sequence $(\mathscr{C}_{n+\nu})$ as formulated
above.

  3. (3)

It follows from the last remark item that the $z$-transform of the sequence
$(J_{n+\nu})$, or the Borel resummation of $\sum_{n}J_{n+\nu}t^{n}$, is a
well-defined analytic functions in two variables as defined in the Definition
2.7.

  4. (4)

We also note that the conclusion of Theorem 3.22 remains valid if the uniform
compacta norm is replaced by any other norm.

######  Theorem 3.24 ((Generating function of Bessel functions)).

Let $\nu\in\mathbb{C}$.

  1. (1)

Let $(\mathscr{C}_{n+\nu})_{n}$ be a bilateral sequence of analytic functions
which is a solution of the Bessel module $\mathcal{B}_{\nu}$ in
$\mathcal{O}^{\mathbb{Z}}$. Then there exists a complex number $C_{\nu}$ such
that

$C_{\nu}\,t^{-\nu}\exp\Big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t}\big{)}\Big{]}\sim\sum_{n=-\infty}^{\infty}\mathscr{C}_{\nu+n}(x)\,t^{n},$
(3.36)

where the symbol “$\sim$” above means that the left-side is the Borel
resummation181818Borel resummation is defined in Definition 2.7. of the right-
hand side.

Moreover, the holonomic system of PDEs (3.14) when realized in
$\mathcal{O}_{dd}$ in Example 2.3 defined by (2.1) is given by

$y_{x}+(1/t-t)/2\,y=0,\quad\nu y+ty_{t}-{x}/2\,(1/t+t)\,y=0.$ (3.37)

To prove Theorem 3.24, we need the following two lemmas, which serve the
purpose of simplifying and solving the holonomic system of PDEs (3.14) in
$\mathcal{O}_{dd}$.

######  Lemma 3.25.

The Bessel module $\mathcal{B}_{\nu}$ is isomorphic to

$\dfrac{\mathcal{A}_{2}}{\mathcal{A}_{2}[\partial_{1}+\frac{1}{2}({1}/{X_{2}}-X_{2})]+\mathcal{A}_{2}[(\nu+X_{2}\partial_{2})-\frac{1}{2}X_{1}(1/X_{2}+X_{2})]}.$
(3.38)

######  Lemma 3.26.

The map

$\begin{array}[]{rcl}\dfrac{\mathcal{A}_{2}}{\mathcal{A}_{2}[\partial_{1}+\frac{1}{2}({1}/{X_{2}}-X_{2})]+\mathcal{A}_{2}[(\nu+X_{2}\partial_{2})-\frac{1}{2}X_{1}(1/X_{2}+X_{2})]}&\stackrel{{\scriptstyle\times
t^{-\nu}\exp[\frac{x}{2}(t-\frac{1}{t})]}}{{\longrightarrow}}&\mathcal{O}_{dd}\end{array}$

is well-defined and left $\mathcal{A}_{2}$-linear. In fact, this linear map is
unique up to a complex multiple.

###### Proof.

It is straightforward to verify that the function
$t^{-\nu}\exp[\frac{x}{2}(t-\frac{1}{t})]$ satisfies the system of PDEs

$\frac{\partial y}{\partial
x}+\frac{1}{2}\big{(}\dfrac{1}{t}-t\big{)}y=0,\qquad\nu y+t\frac{\partial
y}{\partial t}-\frac{x}{2}\big{(}\dfrac{1}{t}+t\big{)}y=0.$ (3.39)

To see the uniqueness, for a general analytic function $y=y(x,\,t)$. The first
equation yields $y(x,t)=g(t)\exp\big{[}(\dfrac{x}{2}(t-\dfrac{1}{t})\big{]}$
for some $g$. Substituting this $y$ into the second differential equation
yields $g(t)=Ct^{-\nu}$ for some $C\in\mathbb{C}$. ∎

Proof of Theorem 3.24. Recall the left $\mathcal{A}_{2}$-linear maps
$\mathfrak{j}$ and $\mathfrak{z}$ as defined in (3.30) and (2.6) respectively.
We have the diagram

${\mathcal{B}_{\nu}}$${\mathcal{O}^{\mathbb{Z}}}$${\widetilde{\mathcal{A}}_{2}}$${\mathcal{O}_{dd}}$$\scriptstyle{\mathfrak{j}}$$\scriptstyle{\times{\mathop{\rm
E}}[\frac{X_{1}}{2}(X_{2}-\frac{1}{X_{2}})]}$$\scriptstyle{\mathfrak{z}}$$\scriptstyle{\times
t^{-\nu}}$ (3.40)

in which
$\widetilde{\mathcal{A}}_{2}:=\mathcal{A}_{2}/\big{[}\mathcal{A}_{2}\partial_{1}+\mathcal{A}_{2}(X_{2}\partial_{2}+\nu)\big{]}$.
Here the left vertical map
$\mathcal{B}_{\nu}\longrightarrow\widetilde{\mathcal{A}}_{2}$ of (3.40) is
given by (3.17) from Theorem 3.16. Since the above diagram commutes up to a
complex multiple, the sum $\sum_{n}\mathscr{C}_{\nu+n}(x)\,t^{n}$, being the
image of $1$ in $\mathcal{O}_{dd}$ via the top-right path, is also a formal
solution to the system of PDEs (3.14). Now the bottom-left path of (3.40)
gives $t^{-\nu}\exp\big{[}\dfrac{x}{2}(t-\dfrac{1}{t})\big{]}$ which is a
solution to (3.14) because of Lemma 3.25 and Lemma 3.26. Since
$\mathcal{B}_{\nu}$ is holonomic, the $\mathbb{C}$-dimension of the local
solution space of (3.14) equals the multiplicity of $\mathcal{B}_{\nu}$, which
is one. So, in summary, we may regard (3.36) holds up to a complex scalar
multiple $C_{\nu}$ which can be identically zero. ∎

The next corollary determines the constant $C_{\nu}$ that appears in (3.36)
for various classes of Bessel functions. Moreover, a connection between these
generating functions and Sonine integrals, see e.g. Watson [64, Chap. VI] is
established.

######  Corollary 3.27.

  1. (1)

Let $\nu\in\mathbb{C}$. Then

$t^{-\nu}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\sim\sum_{n=-\infty}^{\infty}J_{\nu+n}(x)\,t^{n}.$
(3.41)

  2. (2)

Let $\nu\in\mathbb{C}^{\ast}$. Then

$e^{-i\pi/2}\,t^{-\nu}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\sim\sum_{n=-\infty}^{\infty}Y_{\nu+n}(x)\,t^{n}.$
(3.42)

  3. (3)

Let $\nu\in\mathbb{C}$. Then

$t^{-\nu}\exp\big{[}\dfrac{x}{2}\big{(}t+\dfrac{1}{t})\big{]}\sim\sum_{n=-\infty}^{\infty}I_{\nu+n}(x)\,t^{n}.$
(3.43)

  4. (4)

Let $\nu\in\mathbb{C}^{\ast}$. Then

$i\pi\,t^{-\nu}\exp\big{[}-\dfrac{x}{2}\big{(}t+\dfrac{1}{t})\big{]}\sim\sum_{n=-\infty}^{\infty}K_{\nu+n}(x)\,t^{n}.$
(3.44)

where the notation “$\sim$” denotes the Borel resummation as defined in
Definition 2.7.

######  Remark 3.28.

We note that the (3.41) and (3.43) above become well-known equalities when
$\nu$ is an integer.

###### Proof.

  1. (1)

Suppose $\mathscr{C}_{\nu+n}=J_{\nu+n}$ in (3.36). We aim to show that
$C_{\nu}=1$. In particular, we have

$C_{\nu}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}-\frac{J_{0}(x)}{t}\sim\sum_{n=-\infty\atop
n\not=0}^{+\infty}J_{\nu+n}(x)\,t^{n-1}.$ (3.45)

That is, there is a $F\in\mathcal{O}_{dd}$202020 It is clear that
$\sum_{n=-\infty\atop n\not=0}^{+\infty}\frac{J_{\nu+n}(x)}{n}\,t^{n}$ is a
formal primitive of the right-side of (3.45). such that

$C_{\nu}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}-\frac{J_{0}(x)}{t}=\mathfrak{z}\big{(}\partial_{2}(J_{\nu+n}/n)\big{)}=\partial_{2}F$
(3.46)

where the $z$-transform given by the Borel-resummation is defined in (2.6),
and the $\partial_{2}$’s are defined in (2.5), where
$\lim_{t\to-\infty}F(x,t)=0$ if $\mathop{\rm Re}(x)>0$. According to Theorem
2.13, the Borel transform $\mathfrak{B}$ representing the
$\mathfrak{z}$-transform is $\mathcal{A}$-linear, it follows that the left-
side of (3.45) also has a primitive, which can be written as

$\int\Big{\\{}C_{\nu}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}-\frac{J_{\nu}(x)}{t}\Big{\\}}\,dt.$

Let $0<\delta<R$ be given and $|\arg x|<\pi/2$. Consider the contour

$\Gamma_{R,\delta}=(Re^{-i\pi},\,\delta
e^{-i\pi})\cup\mathcal{C}_{\delta}\cup(\delta e^{i\pi},\,Re^{i\pi}),$

where $\mathcal{C}_{\delta}$ denotes the circle centred at the origin with
radius $\delta$. Thus the contour $\Gamma_{R,\delta}$ can be considered a
truncated Hankel’s contour which emanates from $-\infty$ below the negative
real-axis and then back to $-\infty$ above the negative real-axis after
circulating the origin once in an anti-clockwise direction.

Hence

$C_{\nu}\int_{\Gamma_{R,\delta}}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\,dt=\int_{\Gamma_{R,\delta}}\frac{J_{\nu}(x)}{t}\,dt$

We now take the Cauchy principal values on both sides of this equation under
the limits that $R\to+\infty$ and $\delta\to 0$. The right-side yields

$C_{\nu}\int_{\Gamma_{\infty,0}}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\,dt=2\pi
i\,J_{\nu}(x).$

However, since the integral on the left-side converges to a Hankel contour, so
we obtain

$\int_{\Gamma_{R,\delta}}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\,dt\longrightarrow\int_{-\infty}^{(0+)}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\,dt=2\pi
i\,J_{\nu}(x)$

as $R\to+\infty$ and $\delta\to 0$, which actually equals to $2\pi
iJ_{\nu}(x)$ from the classical Schläfli-Sonine integrals which is valid for
$|\arg x|<\frac{\pi}{2}$, see e.g., Watson [64, Chap. VI, p. 176, (2)]. This
proves that $C_{\nu}=1$ after comparing the above calculations. Alternatively,
one can verify directly that the Hankel-type contour as the limit of $R\to
0,\,\delta\to 0$ in (3.46) satisfies the Bessel equation (1.19) and the
following asymptotic

$\begin{split}\int_{-\infty}^{(0+)}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\,dt&\stackrel{{\scriptstyle
t=\frac{2x}{u}}}{{=}}\frac{x^{\nu}}{2^{\nu}}\int_{-\infty}^{(0+)}u^{-\nu-1}\exp\big{[}u-\dfrac{x^{2}}{4u}\big{]}\,du\\\
&\approx\frac{x^{\nu}}{2^{\nu}}\frac{2^{\nu}}{\Gamma(\nu+1)}=\frac{2^{\nu}}{\Gamma(\nu+1)}\end{split}$

when $x\to 0$. See also Watson [64, pp. 175-176]. Hence

$\int_{-\infty}^{(0+)}t^{-\nu-1}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\,dt=2\pi
i\,J_{\nu}(x).$

Therefore, we deduce $C_{\nu}=1$. This proves (3.41).

  2. (2)

Recall that

$Y_{\nu+n}(x)=\frac{\cos(\nu+n)\pi\,J_{n+\nu}-J_{-\nu-n}}{\sin\nu\pi}=\cot\nu\pi\,J_{\nu+n}(x)-(-1)^{n}\csc\nu\pi\,J_{-\nu-n}(x).$

Since $\big{(}J_{-\nu-n}(x)\big{)}_{n}$ corresponds to PDEs (3.14) with $\nu$
replaced by $-\nu$. Hence

$t^{\nu}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\sim\sum_{n=-\infty}^{\infty}J_{-\nu-n}(x)\,t^{n}.$

Thus

$\begin{split}\sum_{n=-\infty}^{\infty}Y_{\nu+n}(x)&\sim\sum_{n=-\infty}^{\infty}\big{[}\cot\nu\pi
J_{\nu+n}(x)-(-1)^{n}\csc\nu\pi J_{-\nu-n}(x)\big{]}\,t^{n}\\\
&=\sum_{n=-\infty}^{\infty}\cot\nu\pi
J_{\nu+n}(x)\,t^{n}-\csc\nu\pi\sum_{n=-\infty}^{\infty}(-1)^{n}J_{-\nu-n}(x)\,t^{n}\\\
&=\cot\nu\pi\sum_{n=-\infty}^{\infty}J_{\nu+n}(x)\,t^{n}-\csc\nu\pi\sum_{n=-\infty}^{\infty}J_{-\nu+n}(x)\,(-1/t)^{n}\\\
&=\cot\nu\pi
t^{-\nu}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}-\csc\nu\pi(-1/t)^{\nu}\exp\big{[}\dfrac{x}{2}\big{(}-\dfrac{1}{t}+t)\big{]}\\\
&=\big{(}\cot\nu\pi-(-1)^{\nu}\csc\nu\pi\big{)}t^{-\nu}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\\\
&=e^{-i\pi/2}t^{-\nu}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\end{split}$

as asserted.

  3. (3)

The proofs of (iii) and (iv) can be derived in a similar way as the proof of
(ii) with the definitions $I_{\nu}(x)=J_{\nu}(ix)$ and
$K_{\nu}(x)=\frac{\pi}{2}(I_{-\nu}(x)-I_{\nu}(x))/\sin\nu\pi$.

∎

By an argument using the Cauchy integral formula, we can now ”extract” the
$n$-th coefficient from the generating function (3.41) above to give a way to
explain why the existence of the following integral representation of Bessel
functions avoids using series expansions as in [66, §17.231], [64], see also
[21, §10.10-11].

######  Corollary 3.29 ((Schläfli-Sonine integrals)).

Let $\nu\in\mathbb{C}$ be given and $|\arg x|<\displaystyle\frac{\pi}{2}$.
Then we have

$J_{\nu+n}(x)=\frac{1}{2\pi
i}\int_{-\infty}^{(0+)}t^{-\nu-1-n}\exp\big{[}\dfrac{x}{2}\big{(}t-\dfrac{1}{t})\big{]}\,dt$
(3.47)

holds for every $n\in\mathbb{Z}$.

###  3.6. Generating function of difference Bessel functions

Let us recall the Example 2.4 that $\mathcal{O}_{\vartriangle}$ represents the
space of analytic functions defined on a domain is a left
$\mathcal{A}_{1}$-module with

$(\partial f)(x)=f(x+1)-f(x),\qquad(Xf)(x)=xf(x-1),\quad
f\in\mathcal{O}_{\vartriangle}.$

We will occasionally use the well-known Newton polynomials written in falling
factorial notations:

$(x)_{0}=1,\quad(x)_{n}=\begin{cases}x(x-1)\cdots(x-n+1),&n\geq 1,\\\
1/[(x+1)\cdots(x+n)]&n\leq-1.\end{cases}$

For non-integer $\nu$, we adopt the notation

$(x)_{\nu}=\frac{\Gamma(x+1)}{\Gamma(x+1-\nu)}.$ (3.48)

######  Example 3.30 ((Difference Exponential function)).

Let $a\in\mathbb{C}$. We have seen in Example 3.3 that the Weyl exponential is
any solution of the $D$-module $\mathcal{A}/\mathcal{A}(\partial-a)$. Now in
particular, the difference exponential function $\exp_{\Delta}(\cdot;a)$ is
defined as a solution of $\mathcal{A}/\mathcal{A}(\partial-a)$ in the
$D$-module $\mathcal{O}_{\Delta}$, i.e.,

$\mathcal{A}/\mathcal{A}(\partial-a)\stackrel{{\scriptstyle\times\exp_{\Delta}(\cdot;a)}}{{\longrightarrow}}\mathcal{O}_{\Delta}.$

To obtain a binomial (or “Newton”) series expansion of
$\exp_{\Delta}(\cdot;a)$, we factorize the above left $\mathcal{A}$-linear map
as

$\mathcal{A}/\mathcal{A}(\partial-a)\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\mathcal{A}/\mathcal{A}\partial\stackrel{{\scriptstyle\times
1}}{{\longrightarrow}}\mathcal{O}_{\Delta},$

in which $S=\sum{\frac{a_{k}}{k!}X^{k}}$ as in (3.1). We can illustrate the
above by the commutative diagram

${\mathcal{A}/\mathcal{A}(\partial-a)}$${\mathcal{A}/\mathcal{A}\partial}$${\mathcal{O}_{\Delta}}$$\scriptstyle{\times
S}$$\scriptstyle{\times\exp_{\Delta}(\cdot;a)}$$\scriptstyle{\times 1}$ (3.49)

As a result, we obtain

$\exp_{\Delta}(x;a)=\sum_{k=0}^{\infty}\frac{a^{k}}{k!}x(x-1)\cdots(x-k+1).$

It is well-known that this binomial series converges absolutely to $(a+1)^{x}$
for any complex $x$ and any $|a|<1$, when an appropriate branch of logarithm
is chosen.

######  Example 3.31 ((Difference trigonometric functions)).

We continue the list of producing difference versions of trigonometric
functions from Example 3.4, where the Weyl sine and the Weyl cosine are
solutions of the $D$-module $\mathcal{A}/\mathcal{A}(\partial^{2}+1)$ in a
$D$-module $N$. Now in particular, we define the difference sine function and
the difference cosine function to be solutions of
$\mathcal{A}/\mathcal{A}(\partial^{2}+1)$ in the $D$-module
$\mathcal{O}_{\Delta}$. To obtain binomial series expansions of these two
functions, we use the factorization

$\mathcal{A}/\mathcal{A}(\partial^{2}+1)\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\mathcal{A}/\mathcal{A}\partial\stackrel{{\scriptstyle\times
1}}{{\longrightarrow}}\mathcal{O}_{\Delta}$ (3.50)

in which we respectively take $S=\sin X$ and $S=\cos X$ as in (3.3) and (3.4).
This factorization can be illustrated by the commutative diagram

${\mathcal{A}/\mathcal{A}(\partial^{2}+1)}$${\mathcal{A}/\mathcal{A}\partial}$${\mathcal{O}_{\Delta}}$$\scriptstyle{\times
S}$$\scriptstyle{\times\mathrm{Sin\,}X/\times\mathrm{Cos\,}X}$$\scriptstyle{\times
1}$ (3.51)

Thus we obtain

$\sin_{\Delta}x:=\frac{1}{2i}[(1+i)^{x}-(1-i)^{x}]=\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{(2k-1)!}(x)_{2k-1}$

and

$\cos_{\Delta}x:=\frac{1}{2}[(1+i)^{x}+(1-i)^{x}]=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k)!}(x)_{2k}.$

These two series converge absolutely if and only if either $\mathop{\rm
Re}(x)>0$ or $x=0$, which are unlike the power series of the usual sine and
cosine functions.

######  Example 3.32 ((Difference Bessel functions)).

Let $\mathcal{A}=\mathcal{A}_{1}$ as before. Recall that the Weyl Bessel of
order $\nu$ (3.6) is any solution of the Bessel operator

$L=(X\partial)^{2}+X^{2}-\nu^{2}$

in some $D$-module $N$, in particular, we define the difference Bessel
function of order $\nu$, denoted by $J^{\Delta}_{\nu}(x)$, to be a solution of
$\mathcal{A}/\mathcal{A}L$ in the $D$-module $\mathcal{O}_{\Delta}$. Recall
the factorization from (3.23)

$\mathcal{A}/\mathcal{A}L\stackrel{{\scriptstyle\times
S}}{{\longrightarrow}}\mathcal{A}/\mathcal{A}(X\partial-\nu)\stackrel{{\scriptstyle\times(x)_{\nu}}}{{\longrightarrow}}\mathcal{O}_{\Delta},$
(3.52)

in which

$S=\mathcal{J}_{\nu}=\sum_{k=0}^{\infty}\frac{(-1)^{k}X^{2k}}{2^{\nu+2k}k!\Gamma(\nu+k+1)}.$

Thus we have obtained the Newton series expansion

$J^{\Delta}_{\nu}(x):=S((x)_{\nu})=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{2^{\nu+2k}k!\Gamma(\nu+k+1)}(x)_{\nu+2k}.$
(3.53)

$J^{\Delta}_{\nu}(x)$ solves the second order difference equation (1.18) or in
its equivalent form

$[(x+2)^{2}-\nu^{2}]y(x+2)-(x+2)(2x+3)y(x+1)+2(x+2)(x+1)y(x)=0,$ (3.54)

which is obtained from the action of the element
$L=(X\partial)^{2}+(X^{2}-\nu^{2})$ in $\mathcal{O}_{\Delta}$, see Example
2.4.

When $x\in\mathbb{Z}$, (3.54) and (3.53) are precisely the difference Bessel
equation and the difference Bessel function $J^{\Delta}_{\nu}(x)$ that appear
to be first studied by Bohner and Cuchta [4]212121Bohner and Cuchta [4] called
$J^{\Delta}_{\nu}(x)$ the “Bessel difference function” and the equation (3.54)
the “Bessel difference equation” instead.. Either by replacing $\nu$ in (3.53)
by $-\nu$ or by considering (3.26), we easily obtain

$J^{\Delta}_{-\nu}(x):=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{2^{-\nu+2k}k!\Gamma(-\nu+k+1)}(x)_{-\nu+2k},$
(3.55)

so that $\\{J^{\Delta}_{\nu},J^{\Delta}_{-\nu}\\}$ are $\mathbb{C}$-linearly
independent solutions of (3.54) in $\mathcal{O}_{\Delta}$, provided that
$2\nu\notin\mathbb{Z}$.

We now provide a statement about an half-plane convergence of the
$J^{\Delta}_{\nu}(x)$.

######  Proposition 3.33.

Let $\nu\in\mathbb{C}$ be arbitrary. Then the difference Bessel function
$J^{\Delta}_{\nu}(x)$ (3.53) converges uniformly in each compact subset of
$\mathop{\rm Re}(x)>0$.

A justification of the above statement will be given in the Appendix §A.1.
However, it was demonstrated in [4, p. 1569] that although the series
$J^{\Delta}_{n}(x)$ diverges for $x<-1$, it converges at $x=-1$. Thus, a
natural question is about the convergence in the strip $\\{x:-1<\mathop{\rm
Re}x<-\frac{1}{2}\\}$.

Similar to Example 3.34, we have

######  Example 3.34.

Let $\nu\in\mathbb{C}$ and consider the bilateral sequence of classical Bessel
functions $(J^{\Delta}_{\nu+n})_{n}$. Then the maps

$\begin{array}[]{rcl}\mathfrak{j}_{\Delta}^{+}:\mathcal{B}_{\nu}&\stackrel{{\scriptstyle\times(J^{\Delta}_{\nu+n})}}{{\longrightarrow}}&\mathcal{O}^{\mathbb{Z}}\end{array}$
(3.56)

and

$\begin{array}[]{rcl}\mathfrak{j}_{\Delta}^{-}:\mathcal{B}_{\nu}&\stackrel{{\scriptstyle\times(J^{\Delta}_{-\nu-n})}}{{\longrightarrow}}&\mathcal{O}^{\mathbb{Z}}\end{array}$
(3.57)

are both left $\mathcal{A}_{2}$-linear. More generally, a general solution of
$\mathcal{B}_{\nu}$ in $\mathcal{O}_{\Delta}^{\mathbb{Z}}$ will be denoted by
$\times\mathrm{(}\mathscr{C}^{\Delta}_{\nu+n}\mathrm{)}$. That is,

$\begin{array}[]{rcl}\mathfrak{j}_{\Delta}:\mathcal{B}_{\nu}&\stackrel{{\scriptstyle\times(\mathscr{C}^{\Delta}_{\nu+n})}}{{\longrightarrow}}&\mathcal{O}_{\Delta}^{\mathbb{Z}}.\end{array}$
(3.58)

The following proposition and corollary are direct consequences of Corollary
3.8 with $N=\mathcal{O}_{\Delta}$ and $\mathscr{C}^{\Delta}_{\nu+n}$.

######  Proposition 3.35 ((Delay-difference formulae)).

Let $\nu\in\mathbb{C}$ and $\mathrm{(}\mathscr{C}^{\Delta}_{\nu+n}\mathrm{)}$
be an arbitrary (bilateral sequence of) solution of
$\mathcal{B}_{\nu}\to\mathcal{O}_{\Delta}^{\mathbb{Z}}$ as defined in (3.58).
Then the $\mathrm{(}\mathscr{C}^{\Delta}_{\nu+n}\mathrm{)}$ satisfies the
following delay-difference formulae

$x\Delta\mathscr{C}^{\Delta}_{\nu}(x-1)+\nu\mathscr{C}^{\Delta}_{\nu}(x)-x\mathscr{C}^{\Delta}_{\nu-1}(x-1)=0$
(3.59)

and

$x\Delta\mathscr{C}^{\Delta}_{\nu}(x-1)-\nu\mathscr{C}^{\Delta}_{\nu}(x)+x\mathscr{C}^{\Delta}_{\nu+1}(x-1)=0.$
(3.60)

for each $n$.

The following three-term recurrence is a direct consequence of Corollary 3.10
with $N=\mathcal{O}_{\Delta}$, or directly from the above two corollary.

######  Proposition 3.36 ((Three-term recurrence)).

For every $\nu\in\mathbb{C}$ and every $x\in\mathbb{C}$, we have

$2\nu\mathscr{C}_{\nu}^{\Delta}(x)-x\mathscr{C}_{\nu-1}^{\Delta}(x-1)-x\mathscr{C}_{\nu+1}^{\Delta}(x-1)=0.$
(3.61)

One deduces the following formulae originally derived by Bohner and and Cuchta
[4, Theorems 5-6, Corollary 7] about $J_{\nu}^{\Delta}$ as consequences of the
above when $\mathscr{C}_{\nu+n}^{\Delta}(x)=J_{\nu+n}^{\Delta}(x)$ for all
$n$.

######  Corollary 3.37 (([4, Theorems 5-6, Corollary 7])).

For every $\nu\in\mathbb{C}$ and every $x\in\mathbb{C}$, we have

$\displaystyle x\Delta J^{\Delta}_{\nu}(x-1)+\nu
J^{\Delta}_{\nu}(x)-xJ^{\Delta}_{\nu-1}(x-1)$ $\displaystyle=0,$ (3.62)
$\displaystyle x\Delta J^{\Delta}_{\nu}(x-1)-\nu
J^{\Delta}_{\nu}(x)+xJ^{\Delta}_{\nu+1}(x-1)$ $\displaystyle=0,$ (3.63)
$\displaystyle 2\nu
J_{\nu}^{\Delta}(x)-xJ_{\nu-1}^{\Delta}(x-1)-xJ_{\nu+1}^{\Delta}(x-1)$
$\displaystyle=0.$ (3.64)

###### Proof.

We provide an alternative proof based on the transmutation formulae from
Proposition 3.6. This more direct approach is similar to the derivation of the
Corollary 3.7. Thus it suffices to prove the (3.62). Since the right-side of
(3.7) annihilates the $CJ^{\Delta}_{\nu}(x)$ for some constant $C$, and the
fact that on the Newton series (3.53) gives

$\begin{split}(\partial+\nu/X)J^{\Delta}_{\nu}(x)&=J_{\nu}^{\Delta}(x)-J_{\nu}^{\Delta}(x-1)+\nu
J_{\nu}^{\Delta}(x+1)/(x+1)\\\
&=\frac{\Gamma(x+1)}{2^{\nu}\Gamma(\nu+1)\Gamma(x+2-\nu)}\Big{(}\frac{2}{\nu}+O\big{(}\frac{1}{x}\big{)}\Big{)}\\\
&=\frac{\Gamma(x+1)}{2^{\nu-1}\Gamma(\nu)\Gamma(x+2-\nu)}\Big{(}1+O\big{(}\frac{1}{x}\big{)}\Big{)}\\\
\end{split}$

as $x\approx 0$. On the other hand, according to the left-side of (3.7), the
$(\partial+\nu/X)J^{\Delta}_{\nu}(x)$ must be equal to
$CJ^{\Delta}_{\nu-1}(x)$ for some constant $C$. But

$J_{\nu-1}^{\Delta}(x)\approx\frac{\Gamma(x+1)}{2^{\nu-1}\Gamma(\nu)\Gamma(x+2-\nu)}$

as $x\approx 0$. One immediately deduces that $C=1$ by comparing the above two
asymptotes. This establishes the (3.62). The verification of (3.63) is
similar. One can derive the (3.64) by subtracting the (3.62) and (3.63). ∎

Now we establish that a difference analogue of the generating function
$e^{\frac{x}{2}(t-\frac{1}{t})}$ of the classical Bessel function for
difference Bessel functions (3.21) is given by the function

$\big{[}\frac{1}{2}\big{(}t-\frac{1}{t}\big{)}+1\big{]}^{x}$

and to show that it serves as a generating function to integer order
difference Bessel function to (3.54) in $\mathcal{O}_{\Delta}$.

######  Theorem 3.38.

Let $|\frac{1}{2}(t-\frac{1}{t})|<1$ and $x\in\mathbb{C}$. Then

$\Big{[}\frac{1}{2}(t-\frac{1}{t})+1\Big{]}^{x}=\sum_{n=-\infty}^{\infty}J^{\vartriangle}_{n}(x)\,t^{n}.$
(3.65)

###### Proof.

Since $|\frac{1}{2}(t-\frac{1}{t})|<1$, so the series

$\Big{[}\frac{1}{2}(t-\frac{1}{t})+1\Big{]}^{x}=\sum_{k=0}^{\infty}\binom{x}{k}\big{[}\frac{1}{2}(t-\frac{1}{t})\big{]}^{k}$

is absolutely convergent. Writing

$\Big{[}\frac{1}{2}(t-\frac{1}{t})\Big{]}^{k}=\sum_{s=0}^{k}\binom{k}{s}(-\frac{1}{2t})^{s}(\frac{t}{2})^{k-s}.$

Then

$\Big{[}\frac{1}{2}(t-\frac{1}{t})+1\Big{]}^{x}=\sum_{k=0}^{\infty}(x)_{k}\sum_{s=0}^{k}\frac{1}{(k-s)!s!}(-\frac{1}{2t})^{s}(\frac{t}{2})^{k-s}.$

It is clear that the power of $t$ goes from $-\infty$ to $\infty$, as $k$ goes
to infinity. Then we exchange the order of summations to arrive at

$\sum_{-\infty}^{\infty}t^{n}\sum_{s=0}^{\infty}(-1)^{s}\frac{(x)_{n+2s}}{s!(n+s)!2^{n+2s}}=\sum_{-\infty}^{\infty}J^{\Delta}_{n}(x)t^{n}.$

∎

The series manipulation technique used in the above proof is in the same
spirit of that for the classical Bessel generating function (3.21) presented
in [64, §2.1]. In order to deal with the generating function when $\nu\not=0$,
we first need to study the asymptotic behaviour of $J^{\Delta}_{\nu+n}$ as
$n\to\infty$. We apply the same technique of Banach algebra to
$(J^{\Delta}_{n+\nu})$ as in the proof of Theorem 3.22 with suitable
modification.

######  Theorem 3.39 ((Asymptotic behaviour)).

The sequence $(J^{\Delta}_{n+\nu})$ is a uniformly 1-Gevrey compacta solution
of $\mathcal{B}_{\nu}$ in $\mathcal{O}_{\Delta}^{\mathbb{Z}}$.

###### Proof.

The proof makes use of the Newton transform $\mathfrak{N}$ which is described
in Theorem 5.1. Fix a compact $K\subset\mathbb{C}$ and for each
$f\in\mathcal{O}_{d}$, one defines a norm by

$|\|f|\|=\|\mathfrak{N}(f)\|_{L^{\infty},K},$

where the right-side is the sup norm on the compact set $K$. Then
$(\mathcal{O}_{d},|\|\cdot|\|)$ becomes a Banach algebra which can be verified
easily.

The proof of Theorem 3.22 under this alternative norm shows that there exists
$C>0$ such that

$|\|J_{\nu+n}|\|\leq Cn!(\dfrac{2}{|\|x|\|})^{n}\quad\mbox{for large }n.$

Consequently, for large $n$,

$\|J^{\Delta}_{\nu+n}\|_{L^{\infty},K}=\|\mathfrak{N}(J_{\nu+n})\|_{L^{\infty},K}=|\|J_{\nu+n}|\|\leq
Cn!(\dfrac{2}{|\|x|\|})^{n}.$

∎

######  Remark 3.40.

Similarly, the sequence $(J^{\Delta}_{\nu-n})$ is also 1-Gevrey. The same is
true if $J^{\Delta}$ is replaced by $\mathscr{C}^{\Delta}$. As a consequence,
the $z$-transform of the sequence $(J^{\Delta}_{\nu+n})$ is a well-defined
analytic function in two variables.

######  Theorem 3.41.

Let $\nu\in\mathbb{C}$.

  1. (1)

Let $(\mathscr{C}^{\Delta}_{n+\nu})_{n}$ be a bilateral sequence of analytic
functions which is a solution of the Bessel module $\mathcal{B}_{\nu}$ in
$\mathcal{O}_{\Delta}^{\mathbb{Z}}$. Then there exists a $1$-periodic function
$C_{\nu}$ such that

$C_{\nu}(x)\,t^{-\nu}\big{[}\dfrac{1}{2}(t-\dfrac{1}{t})+1\big{]}^{x}\sim\sum_{n=-\infty}^{\infty}\mathscr{C}_{\nu+n}^{\Delta}(x)\,t^{n},$
(3.66)

where the symbol $\sim$ means that the left-hand side is the Borel resummation
222222Definition 2.7 of the right-hand side whenever it diverges.

  2. (2)

Moreover, the manifestation of the holonomic system of PDEs (3.14) in
$\mathcal{O}_{\Delta d}$ in Example 2.9 defined by (2.4) is given by

$y(x+1,t)-y(x,t)=\frac{1}{2}(t-\frac{1}{t})y(x,t),\quad\nu
y(x,t)+ty_{t}(x,t)=\frac{x}{2}\big{(}t+\frac{1}{t}\big{)}y(x-1,t).$ (3.67)

To prove the above theorem, we need the following lemma, which serves the
purpose of solving the holonomic system of PDEs that arise from
$\mathcal{B}_{\nu}$ in $\mathcal{O}_{\Delta d}$.

######  Lemma 3.42.

$\begin{array}[]{rcl}\dfrac{\mathcal{A}_{2}}{\mathcal{A}_{2}[\partial_{1}+\frac{1}{2}(1/{X_{2}}-X_{2})]+\mathcal{A}_{2}[2(\nu+X_{2}\partial_{2})-X_{1}(X_{2}+1/{X_{2}})]}&\stackrel{{\scriptstyle\times
t^{-\nu}[\frac{1}{2}(t-\frac{1}{t})+1]^{x}}}{{\longrightarrow}}&\mathcal{O}_{\Delta
d}\end{array}$

is a left $\mathcal{A}_{2}$-linear map.

###### Proof.

It follows from the two generators of the $\mathcal{B}_{\nu}$ defined above
that we solve the system of PDEs

$\partial_{1}y+\frac{1}{2}(\dfrac{1}{X_{2}}-X_{2})y=0,\qquad(\nu+X_{2}\partial_{2})y-\frac{X_{1}}{2}(\dfrac{1}{X_{2}}+X_{2})y=0$

for an analytic function $y=y(x,\,t)$ for the current choices of
$\partial_{j},\,X_{j},\ j=1,\,2$. That is,

$y(x+1,t)-y(x,t)=\frac{1}{2}(t-\frac{1}{t})y(x,t),\qquad\nu
y(x,t)+t\frac{\partial y}{dt}(x,t)=\frac{x}{2}(t+\frac{1}{t})y(x-1,t).$ (3.68)

The first equation yields
$y(x,t)=g(t)\big{[}\dfrac{1}{2}(t-\dfrac{1}{t})+1\big{]}^{x}$ for some $g$.
Substituting this $y$ into the second differential-difference equation yields
$g(t)=C_{\nu}t^{-\nu}$ for some $C_{\nu}\not=0$. ∎

Proof of Theorem 3.41. Recall the maps $\mathfrak{j}_{\Delta}$ and
$\mathfrak{z}_{\Delta}$ as defined in (3.58) and (2.7) respectively. We have
the diagram

${\mathcal{B}_{\nu}}$${\mathcal{O}^{\mathbb{Z}}_{\Delta}}$${\widetilde{A}_{2}}$${\mathcal{O}_{\Delta
d}}$$\scriptstyle{\mathfrak{j}_{\Delta}}$$\scriptstyle{\times{\mathop{\rm
E}}\big{[}\frac{X_{1}}{2}\big{(}X_{2}-\frac{1}{X_{2}}\big{)}\big{]}}$$\scriptstyle{\mathfrak{z}_{\Delta}}$$\scriptstyle{\times
t^{-\nu}}$ (3.69)

in which
$\widetilde{A}_{2}:=\overline{\mathcal{A}_{2}/\big{[}\mathcal{A}_{2}\partial_{1}+\mathcal{A}_{2}(X_{2}\partial_{2}+\nu)\big{]}}$.
Since the above diagram commutes up to a periodic multiple, the sum
$\sum_{n}\mathscr{C}_{\nu+n}^{\Delta}(x)\,t^{n}$, being the image of $1$ in
$\mathcal{O}_{\Delta d}$ via the top right path, is also a formal solution to
the system of PDEs (3.14). Now
$t^{-\nu}\big{[}\dfrac{1}{2}(t-\dfrac{1}{t})+1\big{]}^{x}$ is a solution to
(3.14) because of Lemma 3.42. Since $\mathcal{B}_{\nu}$ is holonomic, the
$\mathbb{C}$-dimension of the local solution space of (3.14) equals the
multiplicity of $\mathcal{B}_{\nu}$, which is one. So (3.66) holds up to a
complex $1$-periodic function $C_{\nu}(x)$. ∎

######  Theorem 3.43.

  1. (1)

Let $\mathop{\rm Re}(x)<\mathop{\rm Re}(\nu)$. Then

$e^{i\pi\nu}\frac{\sin(x-\nu)\pi}{\sin(\pi
x)}\,t^{-\nu}\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sim}\sum_{n=-\infty}^{\infty}J^{\Delta}_{n+\nu}(x)\,t^{n},$
(3.70)

where the notation “$\sim$” means that the left-side is the Borel-resummation
of the right-side.

  2. (2)

Moreover, the the generating function from part (i) satisfies the system of
PDE (3.68)232323See also the Appendix C.5.

###### Proof.

From Theorem 3.41, we have

$C_{\nu}(x)t^{-\nu}\big{[}\dfrac{1}{2}\big{(}t-\dfrac{1}{t})+1\big{]}^{x}\sim\sum_{n=-\infty}^{\infty}J^{\Delta}_{\nu+n}(x)\,t^{n}.$

Subtracting $J_{\nu}(x)$ from both sides and then dividing both sides by $t$,
we obtain

$C_{\nu}(x)t^{-\nu-1}\big{[}\dfrac{1}{2}\big{(}t-\dfrac{1}{t})+1\big{]}^{x}-\frac{J_{\nu}(x)}{t}\sim\sum_{n=-\infty,n\neq
0}^{\infty}J^{\Delta}_{\nu+n}(x)\,t^{n-1}.$ (3.71)

Now note that the right-side of (3.71) has an anti-$\partial_{2}$ (e.g.,
$\partial_{2}(\sum_{n=-\infty,n\neq
0}^{\infty}\frac{J^{\Delta}_{\nu+n}(x)}{n}\,t^{n})$). Since the Borel
resummation is $\mathcal{A}$-linear, it follows that the LHS of (3.71) has an
antiderivative.

Now let $0<\delta<R$ be given and $|\arg x|<\pi/2$, and consider the contour

$\Gamma_{R,\delta}=(Re^{-i\pi},\,\delta
e^{-i\pi})\cup\mathcal{C}_{\delta}\cup(\delta e^{i\pi},\,Re^{i\pi}),$

where $\mathcal{C}_{\delta}$ denotes the circle centred at the origin with
radius $\delta$. Thus the contour $\Gamma_{R,\delta}$ can be considered a
truncated Hankel’s contour which emanates from $-\infty$ below the negative
real-axis and then back to $-\infty$ above the negative real-axis after
circulating the origin once in an anti-clockwise direction. By Cauchy integral
theorem, the integral of the left-side of (3.71) along $\Gamma_{R,\delta}$ is
zero, i.e.,

$C_{\nu}(x)\int_{\Gamma_{R,\delta}}t^{-\nu-1}\big{[}\dfrac{1}{2}\big{(}t-\dfrac{1}{t})+1\big{]}^{x}\,dt=\int_{\Gamma_{R,\delta}}\frac{J^{\Delta}_{\nu}(x)}{t}\,dt=2\pi
iJ^{\Delta}_{\nu}(x).$ (3.72)

By (3.73) which is an analytic result to be obtained in the next subsection,
we have

$C_{\nu}(x)=e^{i\pi\nu}\frac{\sin(x-\nu)\pi}{\sin(\pi x)}.$

∎

We remark that when $\nu=0$ in the above theorem, one can recover the infinite
sum in (3.66) from the generating function
$\big{[}\dfrac{1}{2}(t-\dfrac{1}{t})+1\big{]}^{x}$ by Theorem 3.38, thus
showing that $C_{0}(x)\equiv 1$.

###  3.7. Integral representations

Now we aim to find an integral representation of the difference Bessel
functions that is analogous to the Sonine integral representation in Corollary
3.29.

######  Theorem 3.44.

Let $\nu\in\mathbb{C}$ and $\mathop{\rm Re}(x)<\mathop{\rm Re}(\nu)$. Let $U$
be a Hankel-type contour in the $t$-plane that starts from, $+\infty$ above
the real-axis, circles around the origin in the counter-clockwise direction,
and returns to $+\infty$ below the real-axis, and in particular

$t^{-\nu}\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}\bigg{|}_{\partial
U}=0.$

Then the integral

$y_{\nu}(x)=\int_{U}t^{-\nu-1}\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}\,dt$

  1. (1)

satisfies the delay-difference equations (3.62) and (3.63) above,

  2. (2)

and, in particular, the integral solves the Bessel difference equation of
order $\nu$ (3.54):

$x(x-1)\triangle^{2}y(x-2)+x\triangle y(x-1)+x(x-1)y(x-2)-\nu^{2}y(x)=0.$

###### Proof.

We first note that because of $\mathop{\rm Re}(x)>\mathop{\rm Re}(-\nu)$, so
that integral is a well-defined function of $x$.

  1. (1)

We rewite the (3.62) and (3.63) in the equivalent forms

$\begin{split}&2\triangle y_{\nu}(x)-y_{\nu-1}(x)+y_{\nu+1}(x)=0,\\\ &2\nu
y_{\nu}(x)-xy_{\nu-1}(x-1)-xy_{\nu+1}(x-1)=0.\end{split}$

Then, it is easy to see that

$2\triangle
y_{\nu}(x)=\int_{U}(t^{-\nu}-t^{-\nu-2})\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}dt=y_{\nu-1}(x)-y_{\nu+1}(x)$

holds. To verify the second equation, we substitute the $y_{n}$ into the
equation and apply integration-by-parts once to get

$\begin{split}xy_{\nu-1}(x-1)+xy_{\nu+1}(x-1)&=x\int_{U}(t^{-\nu}+t^{-\nu-2})\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x-1}dt\\\
&=2\int_{U}t^{-\nu}d\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}\\\
&=-2\int_{U}\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}dt^{-\nu}\\\ &=2\nu
y_{\nu}(x),\end{split}$

as required

  2. (2)

The proof to this part actually follows directly from the part (i) and
Proposition 3.14. However, we offer a direct verification. Applying
$\triangle$ to the (3.62) and replacing $\nu$ by $\nu-1$ in the (3.63)
respectively yield

$\begin{split}&x\triangle^{2}y_{\nu}(x-1)+\triangle y_{\nu}(x)-x\triangle
y_{\nu-1}(x-1)-y_{\nu-1}(x)+\nu\triangle y_{\nu}(x)=0,\\\ &x\triangle
y_{\nu-1}(x-1)+xy_{\nu}(x-1)-(\nu-1)y_{\nu-1}(x)=0.\end{split}$

Adding the above two equations yields

$x\triangle^{2}y_{\nu}(x-1)+\triangle y_{\nu}(x)+xy_{\nu}(x-1)-\nu
y_{\nu-1}(x)+\nu\triangle y_{\nu}(x)=0.$

Replacing $x$ by $x-1$, and multiplying $x$ throughout both sides of the above
equation yield

$x(x-1)\triangle^{2}y_{\nu}(x-2)+x\triangle
y_{\nu}(x-1)+x(x-1)y_{\nu}(x-2)-x\nu y_{\nu-1}(x-1)+x\nu\triangle
y_{\nu}(x-1)=0.$

Multiplying by $\nu$ on both sides of (3.62) and substitute the resulting
equation into the last expression yield precisely the Bessel difference
equation of order $\nu$ (3.54).

∎

######  Theorem 3.45 ((Difference Schläfli-Sonine integrals)).

Let $\mathop{\rm Re}(x)<\mathop{\rm Re}(\nu)$. Let $U$ be an Hankel-type
contour parametrised by $t$ that starts at $+\infty$, above the real-axis,
circles around the origin in the counter-clockwise direction, and returns to
$+\infty$ below the real-axis, and such that $|2t+t^{2}|\geq 1$ on $U$. Then

$J^{\Delta}_{\nu}(x)=e^{i\pi\nu}\dfrac{\sin(x-\nu)\pi}{\sin(\pi
x)}\dfrac{1}{2\pi
i}\int_{\infty}^{(0+)}t^{-\nu-1}\Big{[}\frac{1}{2}(t-\frac{1}{t})+1\Big{]}^{x}dt$
(3.73)

###### Proof.

The integrand has branch points $0,\,-1\pm\sqrt{2}$ and $\infty$ when $\nu,x$
are not integers. We pick the branch cuts from $[-1-\sqrt{2},\,0]$ and
$[-1+\sqrt{2},\,+\infty]$ on the real-axis. Since $|2t+t^{2}|\geq 1$ on $U$,
the integrand allows an absolutely convergent expansion

$\begin{split}t^{-\nu-1}\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}&=t^{-\nu-1}(1+\frac{t}{2})^{x}\big{[}1-\frac{1}{2t(1+t/2)}\big{]}^{x}\\\
&=t^{-\nu-1}(1+\frac{t}{2})^{x}\sum_{k=0}^{\infty}\binom{x}{k}(-1)^{k}[\frac{1}{2t(1+t/2)}]^{k}.\\\
\end{split}$

Observe that the lowest term of (3.53) has

$\frac{2^{-\nu}\Gamma(x+1)}{\Gamma(\nu+1)\Gamma(x+1-\nu)}\approx\frac{1}{2^{\nu}\Gamma(\nu+1)\Gamma(1-\nu)},$

as $x\approx 0$.

Next we consider the lowest term of

$\begin{split}&\dfrac{1}{2\pi
i}\int_{\infty}^{(0+)}t^{-\nu-1}\big{[}\frac{1}{2}(t-\frac{1}{t})+1\big{]}^{x}dt\\\
&=\dfrac{1}{2\pi
i}\int_{\infty}^{(0+)}t^{-\nu-1}(1+\frac{t}{2})^{x}\sum_{k=0}^{\infty}\binom{x}{k}(-1)^{k}[\frac{1}{2t(1+t/2)}]^{k}.\\\
\end{split}$

That is,

$\dfrac{1}{2\pi i}\int_{\infty}^{(0+)}t^{-\nu-1}(1+\frac{t}{2})^{x}dt.$

Let $\frac{t}{2}=e^{i\pi}u$ in above integration. Then we have

$\dfrac{1}{2\pi i}e^{-\nu\pi
i}2^{-\nu}\int_{-\infty}^{(0+)}u^{-\nu-1}(1-u)^{x}\,du.$

We recall a classical result about beta function, see for example [63, p. 104
(12)], in an Hankel-type contour which happens to match exactly with the
change of variable $t=2e^{i\pi}u$ above, namely that

$\dfrac{\Gamma(p+q+1)}{\Gamma(p+1)\Gamma(q+1)}=\dfrac{1}{2\pi
i}\int_{-\infty}^{(0+)}t^{-p-1}(1-t)^{-q-1}dt,\quad\mathop{\rm Re}(p+q+1)>0$

holds. Hence we obtain

$\dfrac{1}{2\pi i}e^{-\nu\pi
i}2^{-\nu}\int_{-\infty}^{(0+)}u^{-\nu-1}(1-u)^{x}du=e^{-\nu\pi
i}2^{-\nu}\dfrac{\Gamma(\nu-x)}{\Gamma(\nu+1)\Gamma(-x)}$

which is valid for $\mathop{\rm Re}(x-\nu)<0$. Since

$\dfrac{\Gamma(\nu-x)}{\Gamma(-x)}=\dfrac{\sin(-x)\pi\,\Gamma(x+1)}{\sin(\nu-x)\pi\,\Gamma(x+1-\nu)}.$

Hence we obtain

$e^{i\pi\nu}\dfrac{\sin(x-\nu)\pi}{\sin(\pi x)}\dfrac{1}{2\pi
i}\int_{\infty}^{(0+)}t^{-\nu-1}(1+\frac{t}{2})^{x}dt=\frac{2^{-\nu}\Gamma(x+1)}{\Gamma(\nu+1)\Gamma(x+1-\nu)}\approx\frac{2^{-\nu}}{\Gamma(\nu+1)\Gamma(1-\nu)},$

as $x\approx 0$ which matches the coefficient of the lowest term of the other
side. This implies that both sides of (3.73) are the same and (3.73) holds. ∎

We can modify the Hankel-type contour $U$ in the integral considered above to
obtain a second linearly independent solution to the difference Bessel
equation (3.54).

######  Theorem 3.46.

Let $\mathop{\rm Re}(x)>\mathop{\rm Re}(\nu+1)$. Let $U^{\prime}$ be an
Hankel-type contour parametrised by $s=-1/t$ where $t$ parametrised by the
contour $U$ as is defined in Theorem 3.45. Then

$J^{\Delta}_{-\nu}(x)=e^{i\pi}e^{-2i\pi\nu}\dfrac{\sin(x+\nu)\pi}{\sin\pi
x}\dfrac{1}{2\pi
i}\int_{U^{\prime}}s^{-\nu-1}\Big{[}\frac{1}{2}(s-\frac{1}{s})+1\Big{]}^{x}ds.$
(3.74)

###### Proof.

Notice that both ends of the contour $U^{\prime}$ would be at the origin which
corresponds to $U$ at $+\infty$. The restriction $\mathop{\rm
Re}(x)>\mathop{\rm Re}(\nu+1)$ ensures that the integral converges at the
origin. Substituting $t=-1/s$ into the contour integral yields

$\frac{1}{2\pi
i}\int_{U}(-t)^{\nu-1}\Big{[}\frac{1}{2}(t-\frac{1}{t})+1\Big{]}^{x}dt=\frac{-e^{i\nu\pi}}{2\pi
i}\int_{U}t^{\nu-1}\Big{[}\frac{1}{2}(t-\frac{1}{t})+1\Big{]}^{x}dt$

where we need to be more restrictive under the extra temporary assumption that
$\mathop{\rm Re}(x)<\mathop{\rm Re}(-\nu)$ holds in order to ensure that the
integral converges at $+\infty$. Then the proof of the (3.74) follows closely
that of the last theorem. One then considers the lowest term from the integral
and to compare it with that of the $J^{\Delta}_{-\nu}(x)$ in (3.55). We can
now remove the extra assumption $\mathop{\rm Re}(x)<\mathop{\rm Re}(-\nu)$. ∎

##  4\. Half-Bessel modules I: Bessel polynomials

We recall the formula

$K_{n+\frac{1}{2}}(x)=\sqrt{\frac{\pi}{2}}\,e^{-x}x^{-n-\frac{1}{2}}\theta_{n}(x)$
(4.1)

that connects the modified Bessel function (or Macdonald function) of order
$n+\frac{1}{2}$ and the reversed Bessel polynomials $\theta_{n}(x)$. The
modified Bessel functions $I_{\nu}(x)$ and $K_{\nu}(x)$ are essentially
related to $J_{\nu}(ix)$ and $Y_{\nu}(ix)$ respectively. We refer to Watson
[64] for their precise definitions so that both $I_{\nu}(x)$ and $K_{\nu}(x)$
are real on the real axis. We first modify the adopt the transmutation
formulae for the Bessel modules from subsection §3.2 to half-Bessel modules
that better suit our purpose of discussion in this section. Thus, it is more
natural to discuss the reverse Bessel polynomials $\theta_{n}(x)$ than the
Bessel polynomials $y_{n}(x)=x^{n}\theta_{n}(1/x)$ from our ${D}$-modules
viewpoint. So the discussion of the $y_{n}(x)$ and its difference analogue
will be postponed to the end of this section.

###  4.1. Transmutation formulae

Replacing $\partial$ by $-i\partial$ and $X$ by $iX$, we still have
$[-i\partial,\,iX]=1$. So the transmutation formulae in Proposition 3.6
becomes

######  Lemma 4.1.

For each $\nu\in\mathbb{C}$,

$\displaystyle[(X\partial)^{2}-(X^{2}+(\nu+1)^{2})]\left(\partial-\dfrac{\nu}{X}\right)$
$\displaystyle=\left(\partial-\dfrac{\nu+2}{X}\right)[(X\partial)^{2}-(X^{2}+\nu^{2})],$
(4.2)
$\displaystyle[(X\partial)^{2}-(X^{2}+(\nu-1)^{2})]\left(\partial+\dfrac{\nu}{X}\right)$
$\displaystyle=\left(\partial+\dfrac{\nu-2}{X}\right)[(X\partial)^{2}-(X^{2}+\nu^{2})],$
(4.3)

which induce the left $\mathcal{A}$-linear maps

$\begin{array}[]{rcl}\mathcal{A}/\mathcal{A}((X\partial)^{2}-(X^{2}+(\nu+1)^{2}))&\stackrel{{\scriptstyle\times(\partial-\frac{\nu}{X})}}{{\longrightarrow}}&\mathcal{A}/\mathcal{A}((X\partial)^{2}-(X^{2}+\nu^{2})),\\\
\mathcal{A}/\mathcal{A}((X\partial)^{2}-(X^{2}+(\nu-1)^{2}))&\stackrel{{\scriptstyle\times(\partial+\frac{\nu}{X})}}{{\longrightarrow}}&\mathcal{A}/\mathcal{A}((X\partial)^{2}-(X^{2}+\nu^{2})).\end{array}$

As a result we deduce from this corollary the following well-known formulae
with the help of the asymptotic of $K_{\nu}(x)$ as $x\to 0$.

######  Corollary 4.2.

For each $\nu\in\mathbb{C}$, the modified Bessel functions $K_{\nu}(x)$
satisfies the following recurrence relations.

$\displaystyle K^{\prime}_{\nu}(x)-\dfrac{\nu K_{\nu}(x)}{x}$
$\displaystyle=-K_{\nu+1}(x),$ $\displaystyle K^{\prime}_{\nu}(x)+\dfrac{\nu
K_{\nu}(x)}{x}$ $\displaystyle=-K_{\nu-1}(x).$

The formula (4.1) suggests the following “change of variable” of $X\partial$
to $X\partial-X-\nu$ and keep the $X$ unchanged that the transmutation
formulae in Lemma 4.1 become

######  Proposition 4.3 ((Transmutation formulae)).

For each $\nu\in\mathbb{C}$, we have

$\displaystyle\Big{(}X\partial^{2}-2(\nu+\frac{1}{2}+X)\partial$
$\displaystyle+2(\nu+\frac{1}{2})\Big{)}(X\partial-X-2\nu)$
$\displaystyle=\big{(}X\partial-X-2(\nu+\frac{1}{2})\big{)}\Big{(}X\partial^{2}-2(\nu-\frac{1}{2}+X)\partial+2(\nu-\frac{1}{2})\Big{)},$
(4.4) $\displaystyle\Big{(}X\partial^{2}-2(\nu-\frac{1}{2}+X)\partial$
$\displaystyle+2(\nu-\frac{1}{2})\Big{)}\frac{1}{X}(\partial-1)$
$\displaystyle=\frac{1}{X}(\partial-1-\frac{1}{X})\Big{(}X\partial^{2}-2(\nu+\frac{1}{2}+X)\partial+2(\nu+\frac{1}{2})\Big{)},$
(4.5)

which induce the left $\mathcal{A}$-linear maps

$\begin{array}[]{rcl}\mathcal{A}/\mathcal{A}\Big{(}X\partial^{2}-2(\nu+\frac{1}{2}+X)\partial+2(\nu+\frac{1}{2})\Big{)}&\stackrel{{\scriptstyle\times(X\partial-X-2\nu)}}{{\longrightarrow}}&\mathcal{A}/\mathcal{A}\Big{(}X\partial^{2}-2(\nu-\frac{1}{2}+X)\partial+2(\nu-\frac{1}{2})\Big{)},\\\
\mathcal{A}/\mathcal{A}\Big{(}X\partial^{2}-2(\nu-\frac{1}{2}+X)\partial+2(\nu-\frac{1}{2})\Big{)}&\stackrel{{\scriptstyle\times\frac{1}{X}(\partial-1)}}{{\longrightarrow}}&\mathcal{A}/\mathcal{A}\Big{(}X\partial^{2}-2(\nu+\frac{1}{2}+X)\partial+2(\nu+\frac{1}{2})\Big{)}.\end{array}$

###  4.2. Reverse Bessel polynomial modules

####  4.2.1. Reverse Bessel polynomial operator

Some of Burchnall’s formulae found in [8] are applications of the Weyl-
algebraic formulae to be derived in this section. They are special cases in
the manifestation of $D$-modules $\mathcal{O}_{d}$. We shall change the
manifestation of $D$-modules to $\mathcal{O}_{\Delta}$ when we treat the
difference Bessel polynomials later. We are ready to set $\nu=n+\frac{1}{2}$
and $\nu=n-\frac{1}{2}$ for $n\in\mathbb{Z}$ in (4.4) and (4.5) respectively
to obtain

######  Corollary 4.4.

For each $n\in\mathbb{Z}$, the following are well-defined left
$\mathcal{A}$-linear maps.

$\begin{array}[]{rcl}\mathcal{A}/\mathcal{A}\big{(}X\partial^{2}-2(n+1+X)\partial+2(n+1)\big{)}&\stackrel{{\scriptstyle\times(X\partial-X-2n-1)}}{{\longrightarrow}}&\mathcal{A}/\mathcal{A}\big{(}X\partial^{2}-2(n+X)\partial+2n\big{)},\\\
\mathcal{A}/\mathcal{A}\big{(}X\partial^{2}-2(n-1+X)\partial+2(n-1)\big{)}&\stackrel{{\scriptstyle\times\frac{1}{X}(\partial-1)}}{{\longrightarrow}}&\mathcal{A}/\mathcal{A}\big{(}X\partial^{2}-2(n+X)\partial+2n\big{)}.\end{array}$

Let $\Theta_{0}=1$ and consider for each $n\in\mathbb{Z}$ the gauge
transformations defined by right multiplication by

$G_{n}:=-(X\partial-X-2n-1),\quad n\in\mathbb{Z}.\mod X\partial$

Then we define

$\Theta_{n}(X)=G_{n-1}\,\cdots G_{1}\,G_{0}\,\Theta_{0}(X)\mod X\partial$
(4.6)

to be the reversed Bessel polynomial of degree $n$.

######  Theorem 4.5.

For each $n\in\mathbb{N}\cup\\{0\\}$, the map

$\mathcal{A}/\mathcal{A}\big{(}X\partial^{2}-2(n+X)\partial+2n)\big{)}\stackrel{{\scriptstyle\times\Theta_{n}(X)}}{{\longrightarrow}}\mathcal{A}/\mathcal{A}\big{(}X\partial^{2}-2X\partial\big{)}$

where $\Theta_{n}(X)$ is defined in (4.6) is a well-defined left
$\mathcal{A}$-linear map.

In particular, we have

$\begin{split}\Theta_{1}(X)&=G_{0}\,\Theta_{0}(X)=(-X\partial+X+1)=X+1\quad\mod
X\partial\\\
\Theta_{2}(X)&=G_{1}\,G_{0}\,\Theta_{0}(x)=(-X\partial+X+3)\,(X+1)=X^{2}+3X+3,\quad\mod
X\partial\\\
\Theta_{3}(X)&=G_{2}\,G_{1}\,G_{0}\,\Theta_{0}(x)(-X\partial+x-5)X^{2}+3X+3=X^{3}+6X^{2}+15X+15,\quad\mod
X\partial\\\ \cdots&\cdots\end{split}$ (4.7)

######  Theorem 4.6.

Let $n\in\mathbb{N}_{0}$ and $\Theta_{n}(X)=G_{n-1}\cdots
G_{1}G_{0}\Theta_{0}(X)$. Then we have

$\Theta_{n}(X)=\sum_{k=0}^{n}\frac{(n+k)!}{2^{k}(n-k)!\,k!}X^{n-k},$ (4.8)

and the map

$\mathcal{A}/\mathcal{A}\big{(}X\partial^{2}-2(n+X)\partial+2n)\big{)}\stackrel{{\scriptstyle\times\Theta_{n}(X)}}{{\longrightarrow}}\mathcal{A}/\mathcal{A}\partial$
(4.9)

is left $\mathcal{A}$-linear.

On the other hand, let $\Theta_{0}=1$ and consider the gauge transformation
$H$ defined by a left multiplication by the element
$-\frac{1}{X}(\partial-1)$. Then we have

$1=\Theta_{0}(X)=H^{n}(\Theta_{n}(X))=(H\,\cdots\,H)(\Theta_{n}(X)).$ (4.10)

That is,

$\begin{split}\Theta_{0}(X)&=H(\Theta_{1}(X))=1/X(1-\partial)(X+1)\quad\mod\mathcal{A}\partial\\\
\Theta_{0}(X)&=H^{2}(\Theta_{2}(X))=\big{[}1/X(1-\partial)\big{]}^{2}(X^{2}+3X+3),\quad\mod\mathcal{A}\partial\\\
\Theta_{0}(X)&=H^{3}(\Theta_{3}(X))=\big{[}1/X(1-\partial)\big{]}^{3}(X^{3}+6X^{2}+15X+15),\quad\mod\mathcal{A}\partial\\\
\cdots&\cdots\end{split}$ (4.11)

######  Corollary 4.7.

For each $n$, let $\Theta_{n}$ be defined as above, which is a solution of
$X\partial^{2}-2(n+X)\partial+2n$ in the left $\mathcal{A}$-module
$\mathbb{C}[X]$. Then

$\displaystyle(X\partial-X-2n-1)\Theta_{n}+\Theta_{n+1}=0,$
$\displaystyle(\partial-1)\Theta_{n}+X\Theta_{n-1}=0.$

as elements in $M$.

###### Proof.

We only proof the first expression. The second expression can be dealt with
similarly. Let $\nu=n+\frac{1}{2}$ in (4.4). Then we note that the right-side
of (4.4) annihilates the $\Theta_{N}(X)$, while the left-side of (4.4)
therefore implies that

$(X\partial-X-(2n+1))\Theta_{n}(X)=C\Theta_{n+1}(X),\mod\mathcal{A}\partial$

for some constant $C$. A simple comparison of the highest terms on both sides
yields that $C=-1$. This proves the desired result. ∎

We immediately deduce from the above Corollary the following recursion.

######  Corollary 4.8 ((Three-term recurrence)).

For each $n$, let $\Theta_{n}$ be defined as above, which is a solution of
$X\partial^{2}-2(n+X)\partial+2n$ in the left $\mathcal{A}$-module
$\mathbb{C}[X]$. Then

$\Theta_{n+1}-(2n+1)\Theta_{n}-X^{2}\Theta_{n-1}=0$

as an element in $M$.

We next show that the reverse Bessel polynomial equation [32, p. 8] or more
generally the half-order Bessel polynomial operator

$L_{n}=X\partial^{2}-2(n+X)\partial+2n$

can be studied as a special case of the half-Bessel module below. The
construction of this module is motivated by Corollary 4.7 and the
$\mathcal{A}_{2}$-module structure endowed on $\mathcal{O}^{\mathbb{N}_{0}}$
as in Example 2.15.

######  Definition 4.9.

We define the reverse Bessel polynomial module to be the modified half-Bessel
module to be the left $\mathcal{A}_{2}$-module

$\Theta=\dfrac{\mathcal{A}_{2}}{\mathcal{A}_{2}(\partial_{1}-1+{X_{1}}/{\partial_{2}})+\mathcal{A}_{2}(X_{1}\partial_{1}-2X_{2}\partial_{2}-1-X_{1}+\partial_{2})}.$
(4.12)

We first verify that the Bessel polynomial module $\Theta$ defined above is
holomonic.

######  Theorem 4.10.

The left $\mathcal{A}_{2}$-module $\Theta$ has dimension two and multiplicity
two. In particular, $\Theta$ is holonomic.

###### Proof.

Multiply $\partial_{2}$ throughout the first generator from (4.12) yields

$\partial_{2}\partial_{1}-\partial_{2}+X_{1}$

before eliminating $\partial_{1}$ from the second generator from the (4.12)

$X_{1}\partial_{1}-2X_{2}\partial_{2}-1-X_{1}+\partial_{2}.$

This yields

$\begin{split}\partial_{2}\partial_{1}-\partial_{2}+{X_{1}}&=\partial_{2}\big{(}\frac{2X_{2}}{X_{1}}\partial_{2}-\frac{1}{X_{2}}\partial_{2}+1+1/X_{1}\big{)}-\partial_{2}+X_{1}\\\
&=((2X_{1}-1)\partial_{2}^{2}+3\partial_{2}+X_{1}^{2})/X_{1}\end{split}$
(4.13)

Thus we can replace the original two generators by