text:
# Multiple orthogonal polynomials associated with Macdonald functions

## 1 Macdonald functions

In this paper we will investigate certain polynomials satisfying orthogonality properties with respect to weight functions related to Macdonald functions $`K_\nu (z)`$. These Macdonald functions are modified Bessel functions of the second kind and satisfy the differential equation

$$z^2u^{\prime \prime }+zu^{}(z^2+\nu ^2)u=0,$$

for which they are the solution that remains bounded as $`z`$ tends to infinity on the real line. They can be given by the following integral representations

$`K_\nu (z)`$ $`=`$ $`\left({\displaystyle \frac{\pi }{2z}}\right)^{1/2}{\displaystyle \frac{e^z}{\mathrm{\Gamma }(\nu +1/2)}}{\displaystyle _0^{\mathrm{}}}e^tt^{\nu 1/2}(1+{\displaystyle \frac{t}{2z}})^{\nu 1/2}𝑑t`$ (1.1)
$`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{z}{2}}\right)^\nu {\displaystyle _0^{\mathrm{}}}\mathrm{exp}(t{\displaystyle \frac{z^2}{4t}})t^{\nu 1}𝑑t`$ (1.2)
(see, e.g., The Bateman Manuscript Project \[2, Chapter VII\], Nikiforov and Uvarov \[5, pp. 223–226\]). Useful properties are

$`2K_\nu ^{}(z)`$ $`=`$ $`K_{\nu 1}(z)+K_{\nu +1}(z),`$ (1.3)
$`{\displaystyle \frac{2\nu }{z}}K_\nu (z)`$ $`=`$ $`K_{\nu 1}(z)K_{\nu +1}(z).`$ (1.4)
These functions have the asymptotic behavior

$$K_\nu (z)=\left(\frac{\pi }{2z}\right)^{1/2}e^z[1+O(1/z)],z\mathrm{},$$
(1.5)
and near the origin

$$\{\begin{array}{cc}z^\nu K_\nu (z)=2^{\nu 1}\mathrm{\Gamma }(\nu )+o(1)\hfill & \text{if }z0\hfill \\ K_0(z)=\mathrm{log}z+O(1),\hfill & \text{if }z0\hfill \end{array}$$
(1.6)
The moments of these functions are given explictly by

$$_0^{\mathrm{}}x^\mu K_\nu (x)𝑑x=2^{\mu 1}\mathrm{\Gamma }(\frac{\mu +\nu +1}{2})\mathrm{\Gamma }(\frac{\mu \nu +1}{2}),\mathrm{}(\mu \pm \nu +1)>0.$$
(1.7)
We will use as weight functions the scaled Macdonald functions

$$\rho _\nu (x)=2x^{\nu /2}K_\nu (2\sqrt{x}),x>0.$$
(1.8)
which behave like $`\mathrm{exp}(2\sqrt{x})`$ as $`x\mathrm{}`$. From the properties (1.3) and (1.4) one easily obtains the differential properties

$`(x^\nu \rho _\nu )^{}`$ $`=`$ $`x^{(\nu +1)}\rho _{\nu +1},`$ (1.9)
$`\rho _{\nu +1}^{}`$ $`=`$ $`\rho _\nu ,`$ (1.10)
so that $`\rho _\nu `$ satisfies the second order linear differential equation $`[x^{\nu +1}(x^\nu \rho _\nu )^{}]^{}=\rho _\nu `$. By a simple change of variable in (1.7) we can find the moments of $`\rho _\nu `$ which are given by

$$_0^{\mathrm{}}x^\mu \rho _\nu (x)𝑑x=\mathrm{\Gamma }(\mu +\nu +1)\mathrm{\Gamma }(\mu +1),$$
(1.11)
in particular we see that the $`n`$th moment of $`\rho _0`$ is given by $`(n!)^2`$.

In  A. P. Prudnikov formulated as an open problem the construction of the orthogonal polynomials associated with the weight function $`\rho _\nu `$ (and for more general weights, known as ultra-exponential weight functions). In the present paper we will show that it is more natural to investigate multiple orthogonal polynomials for two Macdonald weights $`\rho _\nu `$ and $`\rho _{\nu +1}`$ since for these multiple orthogonal polynomials one has nice differential properties, a Rodrigues formula, and an explicit recurrence relation.

## 2 Multiple orthogonal polynomials

We will investigate two types of multiple orthogonal polynomials for the system of weights $`(\rho _\nu ,\rho _{\nu +1})`$ $`(\nu 0)`$ on $`[0,\mathrm{})`$ with an additional factor $`x^\alpha `$. Let $`n,m`$ and $`\alpha >1`$, then the type 1 polynomials $`(A_{n,m}^\alpha ,B_{n,m}^\alpha )`$ are such that $`A_{n,m}^\alpha `$ is of degree $`n`$ and $`B_{n,m}^\alpha `$ is of degree $`m`$, and they satisfy the orthogonality conditions

$$_0^{\mathrm{}}[A_{n,m}^\alpha (x)\rho _\nu (x)+B_{n,m}^\alpha (x)\rho _{\nu +1}(x)]x^{k+\alpha }𝑑x=0,k=0,1,2,\mathrm{},n+m.$$
(2.1)
This gives $`n+m+1`$ linear homogeneous equations for the $`n+m+2`$ unknown coefficients of the polynomials $`A_{n,m}^\alpha `$ and $`B_{n,m}^\alpha `$, so that we can find the type 1 polynomials up to a multiplicative factor, which we will fix later. It will be convenient to use the notation

$$q_{n,m}^\alpha (x)=A_{n,m}^\alpha (x)\rho _\nu (x)+B_{n,m}^\alpha (x)\rho _{\nu +1}(x).$$

Type 2 polynomials $`p_{n,m}^\alpha `$ are such that $`p_{n,m}^\alpha `$ is of degree $`n+m`$ and satisfies the multiple orthogonality conditions

$`{\displaystyle _0^{\mathrm{}}}p_{n,m}^\alpha (x)\rho _\nu (x)x^{k+\alpha }𝑑x=0`$ $`k=0,1,2,\mathrm{},n1,`$ (2.2)
$`{\displaystyle _0^{\mathrm{}}}p_{n,m}^\alpha (x)\rho _{\nu +1}(x)x^{k+\alpha }𝑑x=0`$ $`k=0,1,2,\mathrm{},m1.`$ (2.3)
This means that we distribute the $`n+m`$ orthogonality conditions over the two weights $`x^\alpha \rho _\nu `$ and $`x^\alpha \rho _{\nu +1}`$. This gives $`n+m`$ linear homogeneous equations for the $`n+m+1`$ unknown coefficients of $`p_{n,m}^\alpha `$. For type 2 polynomials we will consider monic polynomials, thereby fixing the leading coefficient to be $`1`$ and leaving $`n+m`$ coefficients to be determined from (2.2)–(2.3).

Multiple orthogonal polynomials are related to Hermite-Padé simultaneous rational approximation of a system of Markov functions near infinity, (see, e.g., Nikishin and Sorokin  and Aptekarev ). In the present situation the functions to be approximated are

$$f_1(z)=_0^{\mathrm{}}\frac{x^\alpha \rho _\nu (x)}{zx}𝑑x,f_2(z)=_0^{\mathrm{}}\frac{x^\alpha \rho _{\nu +1}(x)}{zx}𝑑x.$$

Type 1 Hermite-Padé approximation (Latin type) consists of finding polynomials $`A_{n,m}`$, $`B_{n,m}`$, and $`C_{n,m}`$ such that

$$A_{n,m}(z)f_1(z)+B_{n,m}(z)f_2(z)C_{n,m}(z)=O(z^{nm2})z\mathrm{},$$

and type 2 Hermite-Padé approximation (German type) is simultaneous rational approximation of $`(f_1,f_2)`$ with a common denominator $`p_{n,m}`$, i.e., one wants polynomials $`p_{n,m}`$, $`R_{n,m}`$ and $`S_{n,m}`$ such that

$`p_{n,m}(z)f_1(z)R_{n,m}(z)=O(z^{n1}),`$ $`z\mathrm{},`$
$`p_{n,m}(z)f_2(z)S_{n,m}(z)=O(z^{m1}),`$ $`z\mathrm{}.`$
The polynomials $`A_{n,m},B_{n,m}`$ then are precisely the type 1 polynomials $`A_{n,m}^\alpha `$ and $`B_{n,m}^\alpha `$ and the polynomial $`p_{n,m}`$ is the type 2 polynomial $`p_{n,m}^\alpha `$. The numerator $`C_{n,m}`$ is then given by

$$C_{n,m}(z)=_0^{\mathrm{}}\left[\frac{A_{n,m}^\alpha (z)A_{n,m}^\alpha (x)}{zx}\rho _\nu (x)+\frac{B_{n,m}^\alpha (z)B_{n,m}^\alpha (x)}{zx}\rho _{\nu +1}(x)\right]x^\alpha 𝑑x,$$

and for type 2 approximation the numerators are

$`R_{n,m}(z)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{p_{n,m}^\alpha (z)p_{n,m}^\alpha (x)}{zx}}\rho _\nu (x)x^\alpha 𝑑x,`$
$`S_{n,m}(z)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{p_{n,m}^\alpha (z)p_{n,m}^\alpha (x)}{zx}}\rho _{\nu +1}(x)x^\alpha 𝑑x.`$
Observe that we can write the system $`(f_1,f_2)`$ as

$$f_1(z)=_0^{\mathrm{}}f(x)\frac{\rho _{\nu +1}(x)x^\alpha }{zx}𝑑x,f_2(z)=_0^{\mathrm{}}\frac{\rho _{\nu +1}(x)x^\alpha }{zx}𝑑x,$$

where $`f`$ is itself a Markov function

$$f(x)=\frac{\rho _\nu (x)}{\rho _{\nu +1}(x)}=\frac{1}{\pi ^2}_0^{\mathrm{}}\frac{s^1ds}{(x+s)[J_{\nu +1}^2(2\sqrt{s})+Y_{\nu +1}^2(2\sqrt{s})]},$$

which follows from a result by Ismail . Note however that $`f`$ is a Markov function of a positive measure for which not all the moments exist. Nevertheless we can still say that $`(f_1,f_2)`$ is a Nikishin system, which guarantees that the polynomials $`A_{n,m}^\alpha `$, $`B_{n,m}^\alpha `$ and $`p_{n,m}^\alpha `$ all can be computed and their degrees are exactly $`n`$, $`m`$, and $`n+m`$ respectively. Furthermore the zeros of $`q_{n,n}^\alpha `$, $`q_{n+1,n}^\alpha `$, $`p_{n,n}^\alpha `$ and $`p_{n+1,n}^\alpha `$ will all be on $`(0,\mathrm{})`$.

## 3 Differential properties

In this section we will give some differential properties for the type 1 and type 2 multiple orthogonal polynomials. We will only consider the multiple orthogonal polynomials on the diagonal or close to the diagonal, i.e., when $`n=m`$ or $`n=m+1`$. These are the most natural and they satisfy interesting relations.

###### Theorem 1

For the type 2 multiple orthogonal polynomials we have for every $`\alpha >1`$

$$\frac{d}{dx}p_{n,n}^\alpha (x)=2np_{n,n1}^{\alpha +1}(x),\frac{d}{dx}p_{n,n1}^\alpha (x)=(2n1)p_{n1,n1}^{\alpha +1}(x).$$
(3.1)
This means that the differential operator acting on type 2 multiple orthogonal polynomials lowers the degree by one and raises the $`\alpha `$ by one. This should be compared with the corresponding differential property

$$\frac{d}{dx}L_n^\alpha (x)=L_{n1}^{\alpha +1}(x)$$

for Laguerre polynomials (see, e.g., Szegő \[7, Eq. (5.1.14) on p. 102\]). The normalizing constant is different here since Laguerre polynomials are not monic polynomials.

Proof: We begin by using (2.2)

$$_0^{\mathrm{}}p_{n,n}^\alpha (x)\rho _\nu (x)x^{k+\alpha }𝑑x=0,k=0,1,\mathrm{},n1,$$

and then use (1.10) to find

$$_0^{\mathrm{}}p_{n,n}^\alpha (x)\rho _{\nu +1}^{}(x)x^{k+\alpha }𝑑x=0,k=0,1,\mathrm{},n1.$$

Integration by parts then gives

$$p_{n,n}^\alpha (x)\rho _{nu+1}(x)x^{k+\alpha }|_0^{\mathrm{}}+_0^{\mathrm{}}\left[p_{n,n}^\alpha (x)x^{k+\alpha }\right]^{}\rho _{\nu +1}(x)𝑑x=0,k=0,1,\mathrm{},n1.$$

The integrated terms vanish for $`k1`$ and $`\alpha >1`$ because of the behavior of $`\rho _{\nu +1}`$ near $`0`$ and $`\mathrm{}`$, given in (1.5) and (1.6). Now work out the differentiation of the product, then we have

$$\begin{array}{c}_0^{\mathrm{}}[p_{n,n}^\alpha (x)]^{}x^{k+\alpha }\rho _{\nu +1}(x)𝑑x\hfill \\ \hfill +(k+\alpha )_0^{\mathrm{}}p_{n,n}^\alpha (x)x^{k+\alpha 1}\rho _{\nu +1}(x)𝑑x=0,k=1,2,\mathrm{},n1.\end{array}$$

The second integral is zero for $`k=1,2,\mathrm{},n`$ because of the orthogonality condition (2.3), hence we conclude that

$$_0^{\mathrm{}}[p_{n,n}^\alpha (x)]^{}x^{k+\alpha }\rho _{\nu +1}(x)𝑑x=0,k=1,2,\mathrm{},n1,$$

or equivalently

$$_0^{\mathrm{}}[p_{n,n}^\alpha (x)]^{}x^{\mathrm{}+\alpha +1}\rho _{\nu +1}(x)𝑑x=0,\mathrm{}=0,1,\mathrm{},n2.$$
(3.2)
Next we do a similar analysis with (2.3), in which we use (1.9) to find

$$_0^{\mathrm{}}p_{n,n}^\alpha (x)x^{k+\alpha +\nu +1}[x^\nu \rho _\nu (x)]^{}𝑑x=0,k=0,1,\mathrm{},n1.$$

Integration by parts gives

$$\begin{array}{c}p_{n,n}^\alpha (x)x^{k+\alpha +1}\rho _\nu (x)|_0^{\mathrm{}}\hfill \\ \hfill +_0^{\mathrm{}}\left[p_{n,n}^\alpha (x)x^{k+\alpha +\nu +1}\right]^{}x^\nu \rho _\nu (x)𝑑x=0,k=0,1,\mathrm{},n1.\end{array}$$

The integrated terms vanish for $`k0`$ and $`\alpha >1`$. Working out the derivative of the product then gives

$$\begin{array}{c}_0^{\mathrm{}}[p_{n,n}^\alpha (x)]^{}x^{k+\alpha +1}\rho _\nu (x)𝑑x\hfill \\ \hfill +(k+\alpha +\nu +1)_0^{\mathrm{}}p_{n,n}^\alpha (x)x^{k+\alpha }\rho _\nu (x)𝑑x=0,k=0,1,\mathrm{},n1.\end{array}$$

The second integral is zero for $`k=0,1,\mathrm{},n1`$ because of the orthogonality (2.2), hence we have

$$_0^{\mathrm{}}[p_{n,n}^\alpha (x)]^{}x^{k+\alpha +1}\rho _\nu (x)𝑑x=0,k=0,1,\mathrm{},n1.$$
(3.3)
Now $`[p_{n,n}^\alpha ]^{}`$ is a polynomial of degree $`2n1`$ with leading coefficient $`2n`$ (since we normalized the type 2 multiple orthogonal polynomials by taking monic polynomials), and by (3.2) and (3.3) it satisfies the orthogonality conditions (2.2)–(2.3) for the type 2 multiple orthogonal polynomial $`p_{n,n1}^{\alpha +1}`$. By unicity we therefore have $`[p_{n,n}^\alpha (x)]^{}=2np_{n,n1}^{\alpha +1}(x)`$. A similar reasoning also gives the result $`[p_{n,n1}^\alpha (x)]^{}=(2n1)p_{n1,n1}^{\alpha +1}(x)`$. Note however that the analysis does not work when $`m\{n,n1\}`$.

There is a similar differential property for type 1 multiple orthogonal polynomials which complements the differential property of the type 2 multiple orthogonal polynomials given in the previous theorem.

###### Theorem 2

For the type 1 multiple orthogonal polynomials we have for every $`\alpha >0`$

$$\frac{d}{dx}[x^\alpha q_{n,n}^\alpha (x)]=x^{\alpha 1}q_{n+1,n}^{\alpha 1}(x),\frac{d}{dx}[x^\alpha q_{n,n1}^\alpha (x)]=x^{\alpha 1}q_{n,n}^{\alpha 1}(x).$$
(3.4)
This means that the differential operator acting on $`x^\alpha `$ times the type 1 polynomials raises the degree by one and lowers the $`\alpha `$ by one. This should be compared with the corresponding differential property

$$\frac{d}{dx}[e^xx^\alpha L_n^\alpha (x)]=(n+1)e^xx^{\alpha 1}L_{n+1}^{\alpha 1}(x)$$

for Laguerre polynomials. What we will really prove is that the derivative of $`x^\alpha q_{n,n}^\alpha (x)`$ is proportional to $`x^{\alpha 1}q_{n+1,n}^{\alpha 1}(x)`$ (and similarly for the derivative of $`x^\alpha q_{n,n1}^\alpha (x)`$ which is proportional to $`x^{\alpha 1}q_{n,n}^{\alpha 1}(x)`$). The choice of the proportionality factor one in (3.4) will fix the normalization which was left unspecified by the homogeneous system of equations (2.1).

Proof: We begin by using (2.1), which we can write as

$$_0^{\mathrm{}}x^\alpha q_{n,n}^\alpha (x^{k+1})^{}𝑑x=0,k=0,1,\mathrm{},2n.$$

Integration by parts then gives

$$q_{n,n}^\alpha (x)x^{k+\alpha +1}|_0^{\mathrm{}}_0^{\mathrm{}}[x^\alpha q_{n,n}^\alpha (x)]^{}x^{k+1}𝑑x=0,k=0,1,\mathrm{},2n.$$

The integrated terms will vanish for every $`\alpha >1`$ and $`k0`$. For the integrand of the integral we have

$`[x^\alpha q_{n,n}^\alpha (x)]^{}`$ $`=`$ $`[x^\alpha A_{n,n}^\alpha (x)\rho _\nu (x)]^{}+[x^\alpha B_{n,n}^\alpha (x)\rho _{\nu +1}(x)]^{}`$
$`=`$ $`\alpha x^{\alpha 1}[A_{n,n}^\alpha \rho _\nu (x)+B_{n,n}^\alpha (x)\rho _{\nu +1}(x)]`$
$`+x^\alpha [(A_{n,n}^\alpha (x))^{}\rho _\nu (x)+A_{n,n}^\alpha (x)\rho _\nu ^{}(x)`$
$`+(B_{n,n}^\alpha (x))^{}\rho _{\nu +1}(x)+B_{n,n}(x)\rho _{\nu +1}^{}(x)].`$
Now use (1.9)–(1.10) to find

$`[x^\alpha q_{n,n}^\alpha (x)]^{}`$ $`=`$ $`x^{\alpha 1}[((\alpha +\nu )A_{n,n}^\alpha (x)+x(A_{n,n}^\alpha (x))^{}xB_{n,n}^\alpha (x))\rho _\nu (x)`$
$`+(\alpha B_{n,n}^\alpha (x)A_{n,n}^\alpha (x)+x(B_{n,n}^\alpha (x))^{})\rho _{\nu +1}(x)].`$
Observe that $`(\alpha +\nu )A_{n,n}^\alpha (x)+x(A_{n,n}^\alpha (x))^{}xB_{n,n}^\alpha (x)`$ is a polynomial of degree at most $`n+1`$ and $`\alpha B_{n,n}^\alpha (x)A_{n,n}^\alpha (x)+x(B_{n,n}^\alpha (x))^{}`$ is a polynomial of degree at most $`n`$, so that we can write

$$[x^\alpha q_{n,n}^\alpha (x)]^{}=x^{\alpha 1}[P_{n+1}(x)\rho _\nu (x)+Q_n(x)\rho _{\nu +1}(x)]$$
(3.5)
with polynomials $`P_{n+1}`$ and $`Q_n`$ of degree at most $`n+1`$ and $`n`$ respectively. We therefore have

$$_0^{\mathrm{}}x^{\alpha 1}[P_{n+1}(x)\rho _\nu (x)+Q_n(x)\rho _{\nu +1}(x)]x^{k+1}𝑑x=0,k=0,1,\mathrm{},2n.$$

In addition to this, we also have

$$_0^{\mathrm{}}[x^\alpha q_{n,n}^\alpha (x)]^{}𝑑x=x^\alpha q_{n,n}^\alpha (x)|_0^{\mathrm{}}=0,$$

so that

$$_0^{\mathrm{}}[x^\alpha q_{n,n}^\alpha (x)]^{}x^{\mathrm{}}𝑑x=0,\mathrm{}=0,1,\mathrm{},2n+1.$$

In view of (2.1) and (3.5) this means that

$$[x^\alpha q_{n,n}^\alpha (x)]^{}=\mathrm{constant}x^{\alpha 1}q_{n+1,1}^{\alpha 1}(x).$$

Since $`q_{n,n}^\alpha `$ and $`q_{n+1,n}^{\alpha 1}`$ are only determined up to a normalizing factor, we can choose the constant equal to one, thereby fixing the normalization.

We can repeat the analysis for $`q_{n,n1}^\alpha `$ with minor changes. In this case we have

$`[x^\alpha q_{n,n1}^\alpha (x)]^{}`$ $`=`$ $`x^{\alpha 1}[((\alpha +\nu )A_{n,n1}^\alpha (x)+x(A_{n,n1}^\alpha (x))^{}xB_{n,n1}^\alpha (x))\rho _\nu (x)`$
$`+(\alpha B_{n,n1}^\alpha (x)A_{n,n1}^\alpha (x)+x(B_{n,n1}^\alpha (x))^{})\rho _{\nu +1}(x)].`$
and since $`(\alpha +\nu )A_{n,n1}^\alpha (x)+x(A_{n,n1}^\alpha (x))^{}xB_{n,n1}^\alpha (x)`$ is a polynomial of degree at most $`n`$ and $`\alpha B_{n,n1}^\alpha (x)A_{n,n1}^\alpha (x)+x(B_{n,n1}^\alpha (x))^{}`$ a polynomial of degree at most $`n`$ also, the orthogonality

$$_0^{\mathrm{}}[x^\alpha q_{n,n1}^\alpha (x)]^{}x^{\mathrm{}}𝑑x=0,\mathrm{}=0,1,\mathrm{},2n$$

gives the required result.

Observe that the reasoning does not work for $`[x^\alpha q_{n,m}^\alpha (x)]^{}`$ when $`m\{n,n1\}`$.

The proof also shows that

$`A_{n+1,n}^{\alpha 1}(x)`$ $`=`$ $`(\alpha +\nu )A_{n,n}^\alpha (x)+x[A_{n,n}^\alpha (x)]^{}xB_{n,n}^\alpha (x),`$
$`B_{n+1,n}^{\alpha 1}(x)`$ $`=`$ $`\alpha B_{n,n}^\alpha (x)A_{n,n}^\alpha (x)+x[B_{n,n}^\alpha (x)]^{},`$
and similarly

$`A_{n,n}^{\alpha 1}(x)`$ $`=`$ $`(\alpha +\nu )A_{n,n1}^\alpha (x)+x[A_{n,n1}^\alpha (x)]^{}xB_{n,n1}^\alpha (x),`$
$`B_{n,n}^{\alpha 1}(x)`$ $`=`$ $`\alpha B_{n,n1}^\alpha (x)A_{n,n1}^\alpha (x)+x[B_{n,n1}^\alpha (x)]^{}.`$
## 4 Rodrigues formula

As a consequence of Theorem 2 we have the following Rodrigues formula for the type 1 multiple orthogonal polynomials.

###### Theorem 3

The type 1 multiple orthogonal polynomials can be obtained from

$$\frac{d^{2n}}{dx^{2n}}\left(x^{2n+\alpha }\rho _\nu (x)\right)=x^\alpha q_{n,n1}^\alpha (x),\frac{d^{2n+1}}{dx^{2n+1}}\left(x^{2n+1+\alpha }\rho _\nu (x)\right)=x^\alpha q_{n,n}^\alpha (x).$$
(4.1)
Proof: By combining the two formulas in (3.4) we get

$$\frac{d^2}{dx^2}\left(x^\alpha q_{n,n1}^\alpha (x)\right)=x^{\alpha 2}q_{n+1,n}^{\alpha 2}.$$

Iterate this $`k`$ times to get

$$\frac{d^{2k}}{dx^{2k}}\left(x^\alpha q_{n,n1}^\alpha (x)\right)=x^{\alpha 2k}q_{n+k,n+k1}^{\alpha 2k}.$$

Choose $`n=0`$ and $`\alpha =2k+\beta `$ to find

$$\frac{d^{2k}}{dx^{2k}}\left(x^{2k+\beta }q_{0,1}^{2k+\beta }(x)\right)=x^\beta q_{k,k1}^\beta .$$

Now $`q_{0,1}^\alpha (x)=A_{0,1}^\alpha \rho _\nu (x)`$, where $`A_{0,1}^\alpha `$ is a constant, which we will take equal to one to fix normalization. Thus, if we replace $`\beta `$ by $`\alpha `$ and $`k`$ by $`n`$, then we find the first formula in (4.1).

Similarly, we can use

$$\frac{d}{dx}\left(x^\alpha q_{n,n1}^\alpha (x)\right)=x^{\alpha 1}q_{n,n}^{\alpha 1}(x),$$

and differentiate it $`2k`$ times, to find

$$\frac{d^{2k+1}}{dx^{2k+1}}\left(x^\alpha q_{n,n1}^\alpha (x)\right)=x^{\alpha 2k1}q_{n+k,n+k}^{\alpha 2k1}(x).$$

Choose $`n=0`$ and $`\alpha =2k+1+\beta `$ to find

$$\frac{d^{2k+1}}{dx^{2k+1}}\left(x^{2k+1+\beta }\rho _\nu (x)\right)=x^\beta q_{k,k}^\beta (x).$$

Replacing $`\beta `$ by $`\alpha `$ and $`k`$ by $`n`$ then gives the required formula.

The Rodrigues formula allows us to compute the type 1 multiple orthogonal polynomials explicitly. There is a relationship between type 1 and type 2 polynomials (which is not typical for the Macdonald weights but holds in general) and this allows us to compute the type 2 multiple polynomials as well. Indeed, we can write

$$\left(\begin{array}{c}q_{n,n}^\alpha \\ q_{n,n1}^\alpha \end{array}\right)=\left(\begin{array}{cc}A_{n,n}^\alpha & B_{n,n}^\alpha \\ A_{n,n1}^\alpha & B_{n,n1}^\alpha \end{array}\right)\left(\begin{array}{c}\rho _\nu \\ \rho _{\nu +1}\end{array}\right),$$

so that we have

$$\left(\begin{array}{cc}A_{n,n}^\alpha & B_{n,n}^\alpha \\ A_{n,n1}^\alpha & B_{n,n1}^\alpha \end{array}\right)^1\left(\begin{array}{c}q_{n,n}^\alpha \\ q_{n,n1}^\alpha \end{array}\right)=\left(\begin{array}{c}\rho _\nu \\ \rho _{\nu +1}\end{array}\right).$$

Writing the inverse of a matrix as the adjoint matrix divided by the determinant gives

$$\left(\begin{array}{cc}B_{n,n1}^\alpha & B_{n,n}^\alpha \\ A_{n,n1}^\alpha & A_{n,n}^\alpha \end{array}\right)\left(\begin{array}{c}q_{n,n}^\alpha \\ q_{n,n1}^\alpha \end{array}\right)=[A_{n,n}^\alpha B_{n,n1}^\alpha A_{n,n1}^\alpha B_{n,n}^\alpha ]\left(\begin{array}{c}\rho _\nu \\ \rho _{\nu +1}\end{array}\right).$$
(4.2)
The polynomial $`A_{n,n}^\alpha B_{n,n1}^\alpha A_{n,n1}^\alpha B_{n,n}^\alpha `$ is of degree at most $`2n`$ and satisfies

$$\begin{array}{c}_0^{\mathrm{}}[A_{n,n}^\alpha (x)B_{n,n1}^\alpha (x)A_{n,n1}^\alpha (x)B_{n,n}^\alpha (x)]\rho _\nu (x)x^{k+\alpha }𝑑x\hfill \\ \hfill =_0^{\mathrm{}}[B_{n,n1}^\alpha (x)q_{n,n}^\alpha (x)B_{n,n}^\alpha (x)q_{n,n1}^\alpha (x)]x^{k+\alpha }𝑑x=0,k=0,1,\mathrm{},n1,\end{array}$$

where we have used (4.2) and (2.1). Furthermore we also have

$$\begin{array}{c}_0^{\mathrm{}}[A_{n,n}^\alpha (x)B_{n,n1}^\alpha (x)A_{n,n1}^\alpha (x)B_{n,n}^\alpha (x)]\rho _{\nu +1}(x)x^{k+\alpha }𝑑x\hfill \\ \hfill =_0^{\mathrm{}}[A_{n,n1}^\alpha (x)q_{n,n}^\alpha (x)+A_{n,n}^\alpha (x)q_{n,n1}^\alpha (x)]x^{k+\alpha }𝑑x=0,k=0,1,\mathrm{},n1.\end{array}$$

Hence we can conclude that

$$A_{n,n}^\alpha (x)B_{n,n1}^\alpha (x)A_{n,n1}^\alpha (x)B_{n,n}^\alpha (x)=\mathrm{constant}p_{n,n}^\alpha (x),$$

where the constant is the leading coefficient of the polynomial. A similar reasoning using

$$\left(\begin{array}{c}q_{n+1,n}^\alpha \\ q_{n,n}^\alpha \end{array}\right)=\left(\begin{array}{cc}A_{n+1,n}^\alpha & B_{n+1,n}^\alpha \\ A_{n,n}^\alpha & B_{n,n}^\alpha \end{array}\right)\left(\begin{array}{c}\rho _\nu \\ \rho _{\nu +1}\end{array}\right),$$

gives

$$A_{n+1,n}^\alpha (x)B_{n,n}^\alpha (x)A_{n,n}^\alpha (x)B_{n+1,n}^\alpha (x)=\mathrm{constant}p_{n+1,n}^\alpha (x).$$

## 5 Recurrence relation

To simplify the notation we put

$$P_{2n}(x)=p_{n,n}^\alpha (x),P_{2n+1}(x)=p_{n+1,n}^\alpha (x).$$

It is known that the sequence $`P_n(x)`$, $`n=0,1,2,\mathrm{}`$ satisfies a third order recurrence relation of the form

$$xP_n(x)=P_{n+1}(x)+b_nP_n(x)+c_nP_{n1}(x)+d_nP_{n2}(x),$$
(5.1)
(see, e.g.,  or ). Explicit formulas for the recurrence coefficients are given in the following theorem.

###### Theorem 4

The recurrence coefficients in (5.1) are given by

$`b_n`$ $`=`$ $`(n+\alpha +1)(3n+\alpha +2\nu )(\alpha +1)(\nu 1)`$
$`c_n`$ $`=`$ $`n(n+\alpha )(n+\alpha +\nu )(3n+2\alpha +\nu )`$
$`d_n`$ $`=`$ $`n(n1)(n+\alpha 1)(n+\alpha )(n+\alpha +\nu 1)(n+\alpha +\nu ).`$
Proof: We begin with the recurrence coefficients of even index, which are used in the recurrence relation

$$xp_{n,n}^\alpha (x)=p_{n+1,n}^\alpha (x)+b_{2n}p_{n,n}^\alpha (x)+c_{2n}p_{n,n1}^\alpha (a)+d_{2n}p_{n1,n1}^\alpha (x).$$
(5.2)
We will write the polynomials $`p_{n,m}^\alpha `$ explicitly as

$$p_{n,m}^\alpha (x)=\underset{k=0}{\overset{n+m}{}}a_{n,m}^\alpha (k)x^{n+mk},$$

and since we are dealing with monic polynomials, we have

$$a_{n,m}^\alpha (0)=1.$$

Comparing the coefficient of $`x^{2n}`$ in (5.2) gives

$$b_{2n}=a_{n,n}^\alpha (1)a_{n+1,n}^\alpha (1),$$
(5.3)
hence we need to know the coefficients $`a_{n,m}^\alpha (1)`$. Comparing the coefficient of $`x^{2n2}`$ and $`x^{2n3}`$ respectively in (3.1) gives the recurrence

$$(2n1)a_{n,n}^\alpha (1)=2na_{n,n1}^{\alpha +1}(1),(2n2)a_{n,n1}^\alpha (1)=(2n1)a_{n1,n1}^{\alpha +1}(1).$$

Combining these relations gives

$$a_{n,n}^\alpha (1)=\frac{n}{n1}a_{n1,n1}^{\alpha +2}(1),$$

which leads to

$$a_{n,n}^\alpha (1)=na_{1,1}^{\alpha +2n2}(1).$$

In order to obtain an explicit formula, we compute $`p_{1,1}^\alpha `$ explicitly by solving the system of equations

$`{\displaystyle _0^{\mathrm{}}}(x^2+a_{1,1}^\alpha (1)x+a_{1,1}^\alpha (2))x^\alpha \rho _\nu (x)𝑑x`$ $`=`$ $`0`$
$`{\displaystyle _0^{\mathrm{}}}(x^2+a_{1,1}^\alpha (1)x+a_{1,1}^\alpha (2))x^\alpha \rho _{\nu +1}(x)𝑑x`$ $`=`$ $`0.`$
Using the moments (1.11) this gives

$`a_{1,1}^\alpha (1)`$ $`=`$ $`2(2+\alpha )(2+\alpha +\nu ),`$ (5.4)
$`a_{1,1}^\alpha (2)`$ $`=`$ $`(1+\alpha )(2+\alpha )(1+\alpha +\nu )(2+\alpha +\nu ).`$ (5.5)
This gives

$$a_{n,n}^\alpha (1)=2n(\alpha +2n)(\alpha +2n+\nu ),$$
(5.6)
and since $`2na_{n+1,n}^\alpha (1)=(2n+1)a_{n,n}^{\alpha +1}(1)`$ this also gives

$$a_{n+1,n}^\alpha (1)=(2n+1)(\alpha +2n+1)(\alpha +2n+\nu +1).$$
(5.7)
Inserting this in (5.3) gives the requested formula for the recurrence coefficient $`b_{2n}`$.

Next we compare the coefficient of $`x^{2n1}`$ in (5.2) to find

$$c_{2n}=a_{n,n}^\alpha (2)a_{n+1,n}^\alpha (2)b_{2n}a_{n,n}^\alpha (1).$$
(5.8)
This means that we also need to know $`a_{n,n}^\alpha (2)`$ and $`a_{n+1,n}^\alpha (2)`$. Compare the coefficient of $`x^{2n3}`$ and $`x^{2n4}`$ respectively in (3.1), then

$$(2n2)a_{n,n}^\alpha (2)=2na_{n,n1}^{\alpha +1}(2),(2n3)a_{n,n1}^\alpha (2)=(2n1)a_{n1,n1}^{\alpha +1}(2),$$

which combined gives

$$a_{n,n}^\alpha (2)=\frac{(2n)(2n1)}{(2n2)(2n3)}a_{n1,n1}^{\alpha +2}(2)=\frac{(2n)(2n1)}{2}a_{1,1}^{\alpha +2n2}(2).$$

Using (5.5) this gives

$$a_{n,n}^\alpha (2)=n(2n1)(\alpha +2n1)(\alpha +2n)(\alpha +2n+\nu 1)(\alpha +2n+\nu ),$$
(5.9)
and since $`(2n1)a_{n+1,n}^\alpha (2)=(2n+1)a_{n,n}^{\alpha +1}(2)`$ this also gives

$$a_{n+1,n}^\alpha (2)=n(2n+1)(\alpha +2n)(\alpha +2n+1)(\alpha +2n+\nu )(\alpha +2n+\nu +1).$$
(5.10)
Inserting (5.9), (5.10), (5.6) and the expression for $`b_{2n}`$ into (5.8) gives, after some straightforward calculus (or by using Maple) the requested expression for $`c_{2n}`$.

Finally, compare the coefficient of $`x^{2n2}`$ in (5.2), then

$$d_{2n}=a_{n,n}^\alpha (3)a_{n+1,n}^\alpha (3)b_{2n}a_{n,n}^\alpha (2)c_{2n}a_{n,n1}^\alpha (1),$$
(5.11)
so that we need $`a_{n,n}^\alpha (3)`$ and $`a_{n+1,n}^\alpha (3)`$. To this end we compare the coefficient of $`x^{2n4}`$ and $`x^{2n5}`$ respectively in (3.1) to find

$$(2n3)a_{n,n}^\alpha (3)=2na_{n,n1}^\alpha (3),(2n4)a_{n,n1}^\alpha (3)=(2n1)a_{n1,n1}^\alpha (3),$$

which combined gives

$$a_{n+1,n}^\alpha (3)=\frac{2n(2n+1)}{(2n2)(2n3)}a_{n,n1}^{\alpha +2}(3)=\frac{(2n+1)(2n)(2n1)}{6}a_{2,1}^{\alpha +2n2}(3).$$

In order to find $`a_{2,1}^\alpha (3)`$ we will compute $`p_{2,1}^\alpha `$ explicitly by solving the system of equations

$`{\displaystyle _0^{\mathrm{}}}(x^3+a_{2,1}^\alpha (1)x^2+a_{2,1}^\alpha (2)x+a_{2,1}^\alpha (3))x^\alpha \rho _\nu (x)𝑑x`$ $`=`$ $`0`$
$`{\displaystyle _0^{\mathrm{}}}(x^3+a_{2,1}^\alpha (1)x^2+a_{2,1}^\alpha (2)x+a_{2,1}^\alpha (3))x^{\alpha +1}\rho _\nu (x)𝑑x`$ $`=`$ $`0`$
$`{\displaystyle _0^{\mathrm{}}}(x^3+a_{2,1}^\alpha (1)x^2+a_{2,1}^\alpha (2)x+a_{2,1}^\alpha (3))x^\alpha \rho _{\nu +1}(x)𝑑x`$ $`=`$ $`0`$
which by using (1.11) and some calculus gives

$$a_{2,1}^\alpha (3)=(3+\alpha +\nu )(2+\alpha +\nu )(1+\alpha +\nu )(3+\alpha )(2+\alpha )(1+\alpha ).$$
(5.12)
Using this gives

$$\begin{array}{c}a_{n+1,n}^\alpha (3)=\frac{(2n+1)(2n)(2n1)}{6}\hfill \\ \hfill (\alpha +2n+\nu +1)(\alpha +2n+\nu )(\alpha +2n+\nu 1)(\alpha +2n+1)(\alpha +2n)(\alpha +2n1),\end{array}$$
(5.13)
and since $`(2n3)a_{n,n}^\alpha (3)=2na_{n,n1}^{\alpha +1}(3)`$ we also have

$$\begin{array}{c}a_{n,n}^\alpha (3)=\frac{2n(2n1)(2n2)}{6}\hfill \\ \hfill (\alpha +2n+\nu )(\alpha +2n+\nu 1)(\alpha +2n+\nu 2)(\alpha +2n)(\alpha +2n1)(\alpha +2n2).\end{array}$$
(5.14)
Using (5.14), (5.13), (5.9), (5.7) and the formulas for $`b_{2n}`$ and $`c_{2n}`$ in (5.11) then gives the requested expression for $`d_{2n}`$.

In a similar way we proceed with the odd indices which appear in the recurrence relation

$$xp_{n+1,n}^\alpha (x)=p_{n+1,n+1}^\alpha (x)+b_{2n+1}p_{n+1,n}^\alpha (x)+c_{2n+1}p_{n,n}^\alpha (x)+d_{2n+1}p_{n,n1}^\alpha (x).$$
(5.15)
By comparing the coefficients of $`x^{2n+1}`$, $`x^{2n}`$ and $`x^{2n1}`$ we get expressions for $`b_{2n+1}`$, $`c_{2n+1}`$ and $`d_{2n+1}`$ respectively in terms of the coefficients $`a_{n,m}^\alpha `$ and after working out these expressions we get the required formulas.

Observe that the recurrence coefficients have the asymptotic behavior

$$b_n3n^2,c_n3n^4,c_nn^6.$$
(5.16)
```
Department of Mathematics
Katholieke Universiteit Leuven
Celestijnenlaan 200 B
B-3001 Heverlee (Leuven)
Belgium

```


format: text

text:
# Untitled Document

Singular polynomials of

generalized Kasteleyn matrices

Nicolau C. SaldanhaResearch supported by CNPq (Brazil) and ENS-Lyon

Abstract:  Kasteleyn counted the number of domino tilings of a rectangle by considering a mutation of the adjacency matrix: a Kasteleyn matrix $`K`$. In this paper we present a generalization of Kasteleyn matrices and a combinatorial interpretation for the coefficients of the characteristic polynomial of $`KK^{}`$ (which we call the singular polynomial), where $`K`$ is a generalized Kasteleyn matrix for a planar bipartite graph. We also present a $`q`$-version of these ideas and a few results concerning tilings of special regions such as rectangles.

Introduction

Kasteleyn (\[K\]) counted the number of domino tilings of a rectangle by considering a mutation of the adjacency matrix, since then known as a Kasteleyn matrix (\[LL\], \[ST\]). Given a planar bipartite graph $`𝒢`$ there are several Kasteleyn matrices $`K`$ for $`𝒢`$ but, as has been shown independently by David Wilson and Horst Sachs, the singular values of $`K`$ or, equivalently, the eigenvalues of $`KK^{}`$, are independent of the choice of $`K`$. Following a question posed by James Propp (\[P\]), we search for a combinatorial interpretation for these numbers.

In section 1 we introduce generalized Kasteleyn matrices for planar bipartite graphs and present a combinatorial interpretation for the determinant of such matrices $`A`$ in terms of counting matchings. In section 2 we address the main issue of understanding what the coefficients of the characteristic polynomial of $`AA^{}`$ represent, and then, in section 3, we consider the special case of Kasteleyn matrices. In section 4 we present the $`q`$-analogs of these ideas. In section 5 we take a look at rectangles in the plane. Finally, in section 6 we present a few other small examples. We find the language of homology theory helpful and use it throughout the paper. We thank James Propp, Richard Kenyon and Horst Sachs for helpful conversations and emails.

1. Generalized Kasteleyn matrices and their determinants

Let $`𝒢`$ be a planar bipartite graph with $`n`$ white vertices and $`n^{}`$ black vertices. We number the white (black) vertices $`1,2,\mathrm{},n`$ ($`1,2,\mathrm{},n^{}`$). A generalized Kasteleyn matrix for $`𝒢`$ is an $`n\times n^{}`$ complex matrix $`A`$ such that

$$|a_{ij}|=\{\begin{array}{cc}1,\hfill & \text{if the }i\text{-th white vertex and the }j\text{-th black vertex are adjacent,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

Such matrices are conveniently represented by labeling the edges of $`𝒢`$ with complex numbers of norm 1.

We may identify a generalized Kasteleyn matrix $`A`$ with a cocomplex $`𝐀C^1(𝒢,S^1)`$ by making the convention that $`a_{ij}`$ indicates the value of $`𝐀(e_{ij})`$, $`e_{ij}`$ being the oriented edge going from the $`j`$-th black to the $`i`$-th white vertex. A notational confusion must be avoided here: the complex numbers of norm 1 form a multiplicative group but the coefficients for homology or cohomology should be additive groups. Thus, from now on, the symbol $`S^1`$ shall denote the additive group $`R/Z`$ and we denote the exponential $`x\mathrm{exp}(2\pi ix)`$ by $`\eta :S^1C`$. In particular, we write $`a_{ij}=\eta (𝐀(e_{ij}))`$. Since $`C^2(𝒢,S^1)=0`$, any cocomplex $`𝐀`$ is automatically closed and a generalized Kasteleyn matrix $`A`$ defines an element $`𝐚H^1(𝒢,S^1)`$.

There is a natural inclusion $`Z/(2)S^1`$; this defines $`\eta :Z/(2)C`$ with $`\eta (m)=(1)^m`$. We also obtain induced inclusions $`C^1(𝒢,Z/(2))C^1(𝒢,S^1)`$ and $`H^1(𝒢,Z/(2))H^1(𝒢,S^1)`$. For a generalized Kasteleyn matrix $`A`$, $`𝐀C^1(𝒢,Z/(2))`$ if and only if $`A`$ is a real matrix.

Lemma 1.1:  There is a unique element $`𝐤H^1(𝒢,Z/(2))`$ such that for any cycle $`C`$,

$$𝐤(C)m+l+1(mod2)$$
$`(1)`$
where $`m`$ is the number of vertices in the interior of $`C`$ and $`2l`$ is the length of $`C`$.

In the statement above, the word ‘cycle’ is used in the sense of graph theory: $`C`$ is a simple closed curve in the plane composed of edges and vertices of $`𝒢`$ and the interior of $`C`$ is well defined by Jordan’s theorem. Of course, graph theory cycles define homology cycles (i.e., closed elements of $`C_1(𝒢,Z)`$) but the converse is not always true; any homology cycle may nevertheless be written a linear combination of graph theory cycles.

Proof:  Uniqueness is obvious since the above equation gives the value of $`𝐤`$ computed against any cycle and thus, by linearity, against any element of $`H_1(𝒢,Z)`$ (recall that $`H^1(𝒢,Z/(2))=\mathrm{Hom}(H_1(𝒢,Z),Z/(2))`$).

In order to construct $`𝐤`$, we first notice that $`H^1(𝒢,Z/(2))=(Z/(2))^h`$, where $`h`$ is the number of holes (bounded connected components of the complement) of $`𝒢`$: in this identification, the coordinates of $`𝐚`$ corresponding to a given hole is $`𝐚(C)`$, $`C`$ being the outer boundary of the said hole. It is then clear that there exists a unique element of $`H^1(𝒢,Z/(2))`$, which we call $`𝐤`$, satisfying equation $`(1)`$ for all such $`C`$.

Figure 1.1 illustrates a minor complication which has to be kept in mind: the boundary of a hole is not always a cycle in the sense of graph theory. There should be no confusion, however: in (a) we have $`l=2`$, $`m=1`$ for the only hole and in (b) we have $`l=2`$, $`m=0`$ and $`l=2`$, $`m=4`$ for the smaller and bigger hole, respectively.

(a)
(b)
(c)

Figure 1.1

Let $`C`$ be an arbitrary cycle: we prove that equation (1) holds for $`C`$. The interior of $`C`$ minus $`𝒢`$ is a union of holes. If we discard the holes which are completely surrounded by other holes in $`C`$ and consider the outer boundaries $`C_1,\mathrm{},C_k`$ of the remaining ones, we have $`𝐤(C)=_i𝐤(C_i)`$. Since equation (1) holds for each $`C_i`$ we have

$$𝐤(C)k+\underset{i}{}l_i+\underset{i}{}m_i(mod2),$$

where $`l_i`$ and $`m_i`$ correspond to $`C_i`$. Notice that $`L=l+_il_i`$, $`L`$ denoting the number of edges of some $`C_i`$ on $`C`$ or in the interior of $`C`$, and $`M=2l+m_im_i`$, $`M`$ denoting the number of vertices of some $`C_i`$ on $`C`$ or in the interior of $`C`$. Finally, by Euler characteristic, $`kL+M=1`$ and we have equation $`(1)`$ for $`C`$, proving our lemma.

We call $`𝐤H^1(𝒢,Z/(2))`$ as defined in the previous lemma the Kasteleyn class; when the graph $`𝒢`$ is not clear from the context, we write $`𝐤_𝒢`$. The definition of $`𝐤`$ involves $`m`$ and thus appears to depend on the way $`𝒢`$ is drawn in the plane. Indeed, examples (a) and (c) in Figure 1.1 represent equivalent graphs, but the Kasteleyn classes are different; solid lines stand for a label 1 and dashed lines stand for a label $`1`$. A Kasteleyn matrix is a generalized Kasteleyn matrix corresponding to the Kasteleyn class.

We restrict ourselves for the rest of this section to the case $`n=n^{}`$ in order to explore the relationship between matchings and the determinant of square generalized Kasteleyn matrices $`A`$.

It is natural to interpret a matching of $`𝒢`$ as the sum of its edges oriented from black to white and thus as an element of $`C_1(𝒢,Z)`$. The boundary of any matching always equals the sum of all white vertices minus the sum of all black vertices; thus, the difference of two matchings of $`𝒢`$ is closed and may be identified with an element of $`H_1(𝒢,Z)`$. Notice furthermore that the difference of two matchings may be written in a unique way as a sum of disjoint graph theory cycles.

For $`𝐚H^1(𝒢,S^1)`$ and two matchings $`\mu _1`$ and $`\mu _2`$, $`\eta (𝐚(\mu _2\mu _1))`$ is a complex number of absolute value 1. In particular, if $`𝐚H^1(𝒢,Z/(2))`$ then $`\eta (𝐚(\mu _2\mu _1))`$ is 1 or $`1`$. We then say that $`\mu _1`$ and $`\mu _2`$ have the same $`𝐚`$-parity if $`\eta (𝐚(\mu _2\mu _1))=1`$; $`𝐚`$-parity splits the set of matchings into two equivalence classes (occasionally one of these classes may turn out to be empty).

Lemma 1.2:  For any planar graph $`𝒢`$, for any $`𝐚H^1(𝒢,S^1)`$ and for any given matching $`\mu _0`$ we have

$$\underset{\mu _1,\mu _2}{}\eta (𝐚(\mu _2\mu _1))=\left|\underset{\mu }{}\eta (𝐚(\mu \mu _0))\right|^2,$$

where $`\mu _1,\mu _2`$ and $`\mu `$ range over all matchings.

Proof:  We may write the right hand side as

$$\left(\underset{\mu _2}{}\eta (𝐚(\mu _2\mu _0))\right)\left(\underset{\mu _1}{}\overline{\eta (𝐚(\mu _1\mu _0))}\right)=\left(\underset{\mu _2}{}\eta (𝐚(\mu _2\mu _0))\right)\left(\underset{\mu _1}{}\eta (𝐚(\mu _0\mu _1))\right)$$

and distribute to get the left hand side.

Define

$$\delta (𝐚,𝒢)=\underset{\mu _1,\mu _2}{}\eta (𝐚(\mu _1\mu _2)),$$

where $`\mu _1`$ and $`\mu _2`$ range over all matchings; if $`𝒢`$ admits no matchings we define $`\delta (𝐚,𝒢)=0`$. As an example, $`\delta (0,𝒢)`$ is the square of the number of matchings of $`𝒢`$. Also, for $`𝐚H^1(𝒢,Z/(2))`$, $`\delta (𝐚,𝒢)`$ is the square of the difference between the number of matchings in each $`𝐚`$-parity equivalence class.

A matching may also be though of as a bijection from the set of white vertices to the set of black vertices. Thus, if $`\mu _1`$ and $`\mu _2`$ are matchings then $`\mu _1^1\mu _2`$ is a permutation of the set of white vertices: we say that these two matchings have the same permutation parity if and only if the this permutation is even.

Lemma 1.3:  Two matchings have the same permutation parity if and only if they have the same $`𝐤`$-parity, $`𝐤`$ being the Kasteleyn class.

Proof:  Let $`\mu _1`$ and $`\mu _2`$ be two matchings and write $`\mu _1\mu _2`$ as a sum of disjoint cycles $`C_1,\mathrm{},C_N`$ of lengths $`2l_1,\mathrm{},2l_N`$. The interior and exterior of any of these cycles is matchable, thus $`m_1,\mathrm{},m_N`$ as in Lemma 1.1 are all even. From equation $`(1)`$, $`𝐤(C_i)l_i+1(mod2)`$ and thus $`𝐤(\mu _0\mu _1)(l_i+1)(mod2)`$.

The permutation $`\mu _0^1\mu _1`$ can be written as a product of $`N`$ cycles (in the permutation sense) corresponding to $`C_1,\mathrm{},C_N`$ with lengths $`l_1,\mathrm{},l_N`$ and the parity of the permutation $`\mu _1^1\mu _2`$ is thus $`(l_i+1)`$. This proves our claim.

Notice that permutation parity, unlike the Kasteleyn class, does not depend on how $`𝒢`$ is drawn in the plane. A corollary of the previous lemma is thus that if differences of matchings generate $`H_1(𝒢,Z)`$ then the Kasteleyn class of $`𝒢`$ is the same for all planar embeddings.

Lemma 1.4:  For any generalized Kasteleyn matrix $`A`$ we have $`|det(A)|^2=\delta (𝐚+𝐤,𝒢)`$.

Proof:  Each non-zero monomial in the expansion of $`det(A)`$ corresponds to a matching. Thus, each matching $`\mu `$ contributes with a complex number of absolute value 1 to $`det(A)`$. The expression $`\eta (𝐚(\mu \mu _0))`$ obtains, up to a fixed multiplicative constant of absolute value 1, the product of the corresponding elements of $`A`$. From Lemma 1.2, $`𝐤`$-parity is permutation parity, i.e., gives the sign of the monomial in the definition of the determinant. Thus, the contribution of $`\mu `$ to $`det(A)`$ is, again up to a fixed multiplicative constant of absolute value 1, $`\eta ((𝐚+𝐤)(\mu \mu _0))`$, proving our lemma.

As a special case, if $`K`$ is a Kasteleyn matrix, $`|det(K)|`$ is the number of matchings of $`𝒢`$: this is Kasteleyn’s original motivation.

2. Singular polynomials of generalized Kasteleyn matrices

Having provided an interpretation for $`|det(A)|`$ when $`A`$ is square, it is natural to ask about other functions of $`A`$, specially if $`A`$ is not square. We should not expect natural interpretations for the argument of $`det(A)`$ since it depends on the way we assign labels to vertices. Also, a few simple experiments will show that the spectrum of $`A`$ (even if $`A`$ is square) is not a function of $`𝐚`$. The following lemma tells us what functions of $`A`$ are determined by $`𝐚`$.

Lemma 2.1:  Let $`A`$ be a generalized Kasteleyn matrix for $`𝒢`$ and let $`𝐚`$ be the corresponding element of $`H^1(𝒢,S^1)`$. Then the generalized Kasteleyn matrices for $`𝒢`$ also corresponding to $`𝐚`$ are precisely the matrices of the form $`D_1AD_2`$ where $`D_1`$ and $`D_2`$ are unitary diagonal matrices. Furthermore, if $`𝒢`$ is connected, $`D_1AD_2=D_1^{}AD_2^{}`$ if and only if there exists a complex number $`z`$ of absolute value 1 with $`D_1=zD_1^{}`$, $`D_2=z^1D_2^{}`$.

It is possible to give a more elementary proof, but following the spirit of the rest of this paper we phrase the proof in homological language.

Proof:  As we saw in Section 1, generalized Kasteleyn matrices correspond to 1-cocomplexes in $`C^1(𝒢,S^1)`$; two such 1-cocomplexes $`𝐀`$ and $`𝐀^{}`$ induce the same element of $`H^1(𝒢,S^1)`$ if and only if their difference is the coboundary of a 0-cocomplex. A 0-cocomplex $`𝐃`$ is a function assigning an element of $`S^1`$ to each vertex; the $`\eta `$’s of these elements may conveniently be arranged in a pair of unitary diagonal matrices, $`D_w`$ for the white and $`D_b`$ for the black vertices. It is a simple translating process to verify that the cocomplex $`𝐀+d(𝐃)`$ corresponds to the generalized Kasteleyn matrix $`D_wAD_b^1`$, thus proving our first claim. The uniqueness of $`D_1`$ and $`D_2`$ up to a constant multiplicative factor corresponds to the fact that the only closed 0-cocomplexes are the constants, i.e., that $`H^0(𝒢,S^1)=S^1`$ if $`𝒢`$ is connected.

We recall that for any complex $`n\times n^{}`$ matrix $`B`$, there are unitary matrices $`U_1`$ and $`U_2`$ such that $`S=U_1BU_2`$ is a real diagonal matrix with non-increasing non-negative diagonal entries; $`S`$ is well-defined given $`B`$ and its diagonal entries (i.e., the $`s_{ii}`$ entries of $`S`$, even if $`S`$ is not square) are called the singular values of $`B`$. The rows of $`U_1`$ (resp., columns of $`U_2`$) are called the left (resp., right) singular vectors of $`B`$. It is easy to see that the singular values and left (resp. right) singular vectors of $`B`$ are the non-negative square roots of the eigenvalues and the eigenvectors of $`BB^{}`$ (resp., $`B^{}B`$). Inspired in these classical notions, we call the characteristic polynomial of $`BB^{}`$ the singular polynomial of $`B`$: its roots are the squares of the singular values of $`B`$. Also, the singular polynomials of $`B`$ and $`U_1BU_2`$ are equal and the singular polynomials of $`B`$ and $`B^{}`$ differ by a factor of $`t^{nn^{}}`$.

It follows from Lemma 2.1 and the remarks in the previous paragraph that the singular polynomial of $`A`$ is determined by $`𝐚`$: we call it $`P_𝐚`$. The singular values of $`A`$ and, if the singular values are simple, the absolute values of the coordinates of the singular vectors (up to a constant factor) are also determined by $`𝐚`$. We shall now present what we find to be a reasonably natural interpretation for the coefficients of $`P_𝐚`$. While these numbers determine the singular values the question remains whether a nice interpretation exists for the actual singular values and vectors.

Let $`𝒢`$ be a balanced subgraph of $`𝒢`$: the inclusion induces a map $`\pi _{𝒢,}:H^1(𝒢,S^1)H^1(,S^1)`$. More concretely, if $`𝐚`$ corresponds to a generalized Kasteleyn matrix $`A`$ then $`\pi _{𝒢,}(𝐚)`$ corresponds to the submatrix of $`A`$ obtained by picking only the elements for which both row and column correspond to elements of $``$. The simplest interpretation is probably in terms of labels for edges: just keep the old labels. When this causes no confusion, we write $`𝐚`$ instead of $`\pi _{𝒢,}(𝐚)`$: for instance, we write $`\delta (𝐚,)`$ instead of the more correct but cumbersome $`\delta (\pi _{𝒢,}(𝐚),)`$.

Theorem 2.2:  Let $`A`$ be a generalized Kasteleyn matrix and let $`P_𝐚(t)=t^n+a_1t^{n1}+\mathrm{}+a_{n1}t+a_n`$ be the singular polynomial of $`A`$. Then

$$a_m=(1)^m\underset{||=2m}{}\delta (𝐚+𝐤_{},)$$
$`(2)`$
where $``$ ranges over all balanced subgraphs with $`2m`$ vertices.

Notice that for $`m=n`$ equation (2) is equivalent to Lemma 1.4. For $`m=1`$, we get the simple remark that $`|a_1|`$ is the number of edges of $`𝒢`$, regardless of $`𝐚`$. An interpretation for $`a_2`$ is already subtler: each subgraph with two white and two black vertices contributes with a real number between 0 and 4. Subgraphs which are not matchable of course contribute with 0 and two disjoint edges as well as four points on a line contribute with 1. The interesting part are the squares, which admit two matchings, say $`\mu _1`$ and $`\mu _2`$: then $`\delta (𝐚+𝐤_{},)=|1\eta (𝐚(\mu _1\mu _2))|^2`$; in general, this may be any number between 0 and 2 but if $`𝐚H^1(𝒢,Z/(2))`$ then this is 0 or 4.

In order to prove this Theorem, we need an auxiliary result in linear algebra. The proof of Lemma 2.3 (actually a rather straightforward computation) may be found in \[??\].

Lemma 2.3:  Let $`P(t)=t^n+a_1t^{n1}+\mathrm{}+a_{n1}t+a_n`$ be the singular polynomial of $`A`$ (where $`A`$ is an arbitrary $`n\times n^{}`$ complex matrix). Then

$$a_m=(1)^m\underset{B}{}|detB|^2$$

where $`B`$ ranges over all $`m\times m`$ submatrices of $`A`$.

Proof of Theorem 2.2:  Since balanced subgraphs of $`𝒢`$ with $`2m`$ elements correspond to $`m\times m`$ submatrices of $`A`$, this is a consequence of Lemma 1.4 and Lemma 2.3.

We now describe another, more graphical, interpretation for Theorem 2.2. We define a pipe system of $`𝒢`$ as an oriented pair $`\nu =(\mu _1,\mu _2)`$ of matchings of a subgraph $``$ of $`𝒢`$; we call $``$ the support of the pipe system. Figure 2.1 shows an example of a pipe system: we draw the edges of $`\mu _2`$ oriented from black to white and the edges of $`\mu _1`$ from white to black (unused edges are represented by dotted lines). A pipe system is thus a collection of pipes (i.e., oriented edges of $`𝒢`$) such that, at each vertex, there is either one pipe coming in and one pipe going out or no pipe coming in and no pipe going out (the water that comes in must go out and you can not pipe too much water through a vertex). We define the size $`|\nu |`$ of the pipe system as half the number of vertices in $``$ (in Figure 3.3, $`|\nu |=5`$).

The Kasteleyn class $`𝐤_{}`$ shall be called $`𝐤_\nu `$. A pipe system obtains an element of $`C_1(,Z)`$ (and thus of $`C_1(𝒢,Z)`$) but must not be confused with it: if two pipes cancel each other homologically, they still have to be taken into account for the pipe system. If $`𝐚H^1(,S^1)`$, $`\eta (𝐚(\nu ))`$ and $`\eta ((𝐚+𝐤_\nu )(\nu ))`$ are well defined complex numbers.

Figure 2.2

Corollary 2.4:  Let $`A`$ be a generalized Kasteleyn matrix and let $`P_𝐚(t)=t^n+a_1t^{n1}+\mathrm{}+a_{n1}t+a_n`$ be the singular polynomial of $`A`$. Then

$$P_𝐚(t)=\underset{\nu }{}(1)^{|\nu |}t^{n|\nu |}\eta ((𝐚+𝐤_\nu )(\nu )),$$

where $`\nu `$ ranges over all pipe systems of $`𝒢`$.

Proof:  This follows directly from Theorem 2.2 and the definitions.

3. Singular polynomials of planar graphs

Theorem 2.2 provides an interpretation for the coefficients of singular polynomials of arbitrary generalized Kasteleyn matrices. In this section we take a closer look at the right hand side of equation (2) when $`A`$ is a Kasteleyn matrix.

For planar balanced bipartite graphs $`𝒢`$, define $`𝐩_{𝒢,}H^1(,Z/(2))`$ by $`𝐩_{𝒢,}=\pi _{𝒢,}(𝐤_𝒢)𝐤_{}`$. In Figure 3.1 we illustrate the several objects involved in this definition: (a), (b), (c) and (d) represent $`𝐤_𝒢`$, $`\pi _{𝒢,}(𝐤_𝒢)`$, $`𝐤_{}`$ and $`𝐩_{𝒢,}`$, respectively, where again solid lines stand for a label 1 and dashed lines stand for a label $`1`$. The following lemma provides an alternate definition for this class.

Lemma 3.1:  Let $`𝒢`$ be balanced planar graphs and let $`C`$ be a cycle in $``$. Let $`q`$ be the number of vertices of $`𝒢`$ not belonging to $``$ which are inside $`C`$. Then

$$𝐩_{𝒢,}(C)q(mod2).$$

Proof:  This follows directly from equation $`(1)`$ in Lemma 1.1.

(a)
(b)
(c)
(d)

Figure 3.1

In the hope of making the intuitive meaning of this definition clearer, especially for adjacency graphs of quadriculated or triangulated disks, we introduce some extra structure.

Let $`\overline{𝒢}`$ be the CW-complex obtained from $`𝒢`$ by closing each hole with a 2-cell; $`\overline{𝒢}`$ is thus always homeomorphic to a disk. For $`𝒢`$, let $`\overline{}\overline{𝒢}`$ be the obtained from $``$ by adding the 2-cells of $`\overline{𝒢}`$ whose boundaries are contained in $``$; in other words, we close the holes of $``$ which contain no points of $`𝒢`$. The inclusion $`\overline{}`$ induces an injective map from $`H^1(\overline{},Z/(2))`$ to $`H^1(,Z/(2))`$ which allows for a natural identification of $`H^1(\overline{},Z/(2))`$ with a subset of $`H^1(,Z/(2))`$.

Lemma 3.2:  $`𝐩_{𝒢,}`$ belongs to $`H^1(\overline{},Z/(2))`$.

Proof:  This follows easily from Lemma 3.1.

Recall that two tilings by dominoes of a quadriculated region are said to differ by a flip if they coincide except for two dominoes; in other words, their difference (in the homological sense) is a square. If $`𝒢`$ is the graph of a quadriculated planar region, the difference between two tilings of $``$ differing by a flip is 0 in $`H_1(\overline{},Z)`$; thus, tilings mutually accessible by flips always have the same $`𝐩_{𝒢,}`$-parity.

We may now state the promised interpretation for the coefficients of singular polynomial $`P_𝐤`$ of $`K`$. Since $`P_𝐤`$ is well defined from $`𝒢`$, we may adopt a lighter notation and call it $`P_𝒢`$, the singular polynomial of $`𝒢`$.

Theorem 3.3:  Let $`𝒢`$ be a planar bipartite graph and let $`P_𝒢=t^n+k_1t^{n1}+\mathrm{}+k_{n1}t+k_n`$ be the singular polynomial of $`𝒢`$. Then

$$k_m=(1)^m\underset{||=2m}{}\delta (𝐩_{𝒢,},)$$
$`(3)`$
where $``$ ranges over all balanced subgraphs with $`2m`$ vertices.

Proof:  This is a corollary of Theorem 2.2 and the definition of $`𝐩_{𝒢,}`$.

Recall that if $`𝐩_{𝒢,}=0`$ then $`\delta (𝐩_{𝒢,},)`$ is just the square of the number of matchings of $``$. This always happens if $`\overline{}`$ is simply connected. As a corollary, if $`𝒢`$ is the graph of a quadriculated planar region and $`m3`$, or if $`𝒢`$ is the graph of a triangulated planar region and $`m5`$, then

$$k_m=(1)^m\underset{||=2m}{}\delta (0,)$$

where $``$ ranges over all balanced subgraphs (subregions) with $`2m`$ vertices (squares, triangles).

Notice that $`𝐩_{𝒢,}`$, and thus the right hand side of equation (3), depends on the way $`𝒢`$ is drawn in the plane. Examples (a) and (c) in Figure 1.1 show that $`P_𝒢`$ indeed depends on the way $`𝒢`$ is drawn: for (a) we have $`k_2=9`$ but for (b) we have $`k_2=5`$. This causes the singular values to change in a complicated way: for (a) the singular values are approximately $`0.5549581321,0.8019377358,2.246979604`$ while for (c) they are $`0.3472963553,1.532088886,1.879385242`$.

It is natural to conjecture that the number of non-zero singular values coincides with the size of a maximal partial matching of $`𝒢`$. In Figure 3.2(a) we present an example to show that this is not always true: there are partial matchings of size 3 but since the singular polynomial of $`𝒢`$ is $`t^37t^2+10t`$ there are only two non-zero singular values: $`\sqrt{2}`$ and $`\sqrt{5}`$. In Figure 3.2(b) we draw the same graph in a different way and we now have three non-zero singular values: 1, $`\sqrt{2}`$ and 2.

(a)
(b)

Figure 3.2

We state Theorem 3.3 in the language of pipe systems. Denote $`𝐩_{𝒢,}(\nu )`$ (where $``$ is the support of $`\nu `$) by $`𝐩(\nu )Z/(2)`$. We describe an elementary definition of $`𝐩(\nu )`$. Join pairs of vertices not in $``$, matching black vertices with white vertices in an arbitrary way; if all intersections are transversal, $`𝐩(\nu )`$ is the parity of the number of intersections of such new lines with the pipes.

Corollary 3.4:  Let $`𝒢`$ be a planar graph and $`P_𝒢(t)`$ its singular polynomial. Then

$$P_𝒢(t)=\underset{\nu }{}(1)^{|\nu |}t^{n|\nu |}\eta (𝐩(\nu )),$$

where $`\nu `$ ranges over all pipe systems of $`𝒢`$.

Proof:  This follows directly from Theorem 3.3 and Corollary 2.4.

4. $`q`$

In several branches of combinatorics, $`q`$-analogues or quantizations of classical problems have been seen to be interesting and useful. There are often several interpretations for the $`q`$-analogue of a given concept, some sophisticated (involving quantum groups and the like) and some elementary. In this section we briefly consider a $`q`$-analogue of Kasteleyn matrices in a very naïve way and extend the results of the previous sections to this setting; our interest in doing so is that the methods of the previous sections extend very easily to this more general context and the coefficients will actually have a natural interpretation.

Let $`C_q=C[q,q^1]`$. We extend the usual complex conjugation to $`C_q`$ by postulating $`\overline{q}=q^1`$; $`q`$ may be thought of as an unknown complex number of absolute value 1. Let $`𝒢`$ be a planar bipartite graph with $`n`$ white vertices and $`n^{}`$ black vertices. A generalized Kasteleyn $`q`$-matrix for $`𝒢`$ is an $`n\times n^{}`$ matrix $`A`$ with coefficients in $`C_q`$ such that

$$a_{ij}\overline{a}_{ij}=\{\begin{array}{cc}1,\hfill & \text{if the }i\text{-th white vertex and the }j\text{-th black vertex are adjacent,}\hfill \\ 0\hfill & \text{otherwise;}\hfill \end{array}$$

thus, the entries are always monomials and substituting $`q`$ by a complex number of absolute value 1 changes a generalized Kasteleyn $`q`$-matrix into an ordinary generalized Kasteleyn matrix.

Consider the additive group $`S_q=S^1Z`$ and let $`𝐪`$ be the canonical generator of the $`Z`$ component. If we extend the classical $`\eta `$ to $`\eta :S^1ZC_q`$ by postulating $`\eta (𝐪)=q`$ we may identify a generalized Kasteleyn $`q`$-matrix $`A`$ with a cocomplex $`𝐀C^1(𝒢,S_q)`$. Again, $`C^2(𝒢,S_q)=0`$ and $`A`$ defines an element $`𝐚H^1(𝒢,S_q)`$.

Let $`𝒢`$ be a planar graph. Bounded connected components of the complement of $`𝒢`$ have well defined positively oriented boundaries $`\beta `$ in $`H_1(𝒢,Z)`$. We define a Kasteleyn $`q`$-matrix to be a generalized Kasteleyn $`q`$-matrix $`A`$ such that $`𝐚(\beta )=𝐪`$ for all such boundaries $`\beta `$. We define the singular $`q`$-polynomial of $`𝒢`$ to be the singular polynomial of a Kasteleyn $`q`$-matrix of $`𝒢`$. As before, singular $`q`$-polynomials are easuly seen to be well defined but now they are of course polynomials in $`C_q[t]`$, or, rather equivalently, polynomials in two variables $`q`$ and $`t`$. Finally, define the area of a pipe system $`\nu `$, $`A(\nu )`$ to be the number of bounded connected components of the complement of $`𝒢`$ positively surrounded by $`\nu `$, counted with sign and multiplicity.

With these definitions we have the following theorem.

Theorem 4.1:  Let $`𝒢`$ be a planar graph and $`P_𝒢(q,t)`$ its singular $`q`$-polynomial. Then

$$P_𝒢(q,t)=\underset{\nu }{}(1)^{|\nu |}q^{A(\nu )}t^{n|\nu |}\eta (𝐩(\nu )),$$

where $`\nu `$ ranges over all pipe systems of $`𝒢`$.

Since the proof is entirely analogous to that of Corollary 3.4, we leave the details to the reader.

5. Rectangles

Kasteleyn (\[K\]) computes the determinant of $`K`$ for rectangles essentially by computing its singular values. For the reader’s convenience, we repeat that part of his work in our language. In order to simplify notation in the statement and proof, let

$$\begin{array}{cc}\hfill X_{M,N}^+& =\left\{(k,\mathrm{})Z^2\right|1k\frac{M+1}{2}\text{ and if }k=\frac{M+1}{2}\text{ then }1\mathrm{}\frac{N+1}{2}\},\hfill \\ \hfill X_{M,N}^{}& =\left\{(k,\mathrm{})Z^2\right|1k\frac{M+1}{2}\text{ and if }k=\frac{M+1}{2}\text{ then }1\mathrm{}<\frac{N+1}{2}\}.\hfill \end{array}$$

Theorem 5.1:  Let $`𝒢`$ be a $`M\times N`$ rectangular grid and let $`K`$ be its Kasteleyn matrix. Then the non-zero singular values of $`K`$ are $`\sigma _{k,\mathrm{}}`$, $`(k,\mathrm{})X_{M,N}^{}`$, where

$$\sigma _{k,\mathrm{}}^2=(\alpha ^k+\alpha ^k)^2+(\beta ^{\mathrm{}}+\beta ^{\mathrm{}})^2,\alpha =\mathrm{exp}\left(\frac{\pi i}{M+1}\right),\beta =\mathrm{exp}\left(\frac{\pi i}{N+1}\right).$$

The complicated description of the allowed values of the indices $`k`$ and $`\mathrm{}`$ is necessary in order to avoid zeroes and duplications in a way which is correct for all possible parities of $`M`$ and $`N`$ (Kasteleyn has a simpler formula since he assumes $`N`$ to be even). Notice that

$$\sigma _{k,\mathrm{}}=2\left(\mathrm{cos}^2\frac{k\pi }{M+1}+\mathrm{cos}^2\frac{\mathrm{}\pi }{N+1}\right)^{1/2},$$

(an expression closer to Kasteleyn’s), $`\sigma _{k,\mathrm{}}=\sigma _{M+1k,\mathrm{}}=\sigma _{k,N+1\mathrm{}}=\sigma _{M+1k,N+1\mathrm{}}`$ and that $`\sigma _{k,\mathrm{}}=0`$ if and only if $`M`$ and $`N`$ are both odd, $`k=(M+1)/2`$ and $`\mathrm{}=(N+1)/2`$. Thus, if we just demand $`1kM`$ and $`1\mathrm{}N`$ then all non-zero singular values are counted twice and we occasionally introduce a 0.

Proof:  We index vertices by pairs $`(k^{},\mathrm{}^{})`$, $`1k^{}m`$, $`1\mathrm{}^{}n`$. The vertex $`(k^{},\mathrm{}^{})`$ is called white when $`k^{}+\mathrm{}^{}`$ is even. Define $`K`$ as the Kasteleyn matrix with entries 1 for horizontal edges and $`i`$ for vertical edges: $`K`$ defines a linear transformation from the “black space” to the “white space”. Consider the white vectors

$$w_{k,l}=(\alpha ^{kk^{}}\alpha ^{kk^{}})(\beta ^{\mathrm{}\mathrm{}^{}}\beta ^{\mathrm{}\mathrm{}^{}}),$$

$`(k,\mathrm{})X_{m,n}^+`$: they clearly form an orthogonal basis for the white space (this is where a careful choice of $`X_{m,n}^+`$ becomes necessary). Similarly, the black vectors $`b_{k,l}`$ defined by the same formula with $`(k,\mathrm{})X_{m,n}^{}`$ form an orthogonal basis for the black space. A simple computation yields $`|w_{k,l}|=|b_{k,l}|`$ for $`(k,l)X_{m,n}^{}`$ and

$$\begin{array}{cc}\hfill Kb_{k,\mathrm{}}& =\left((\alpha ^k+\alpha ^k)+i(\beta ^{\mathrm{}}+\beta ^{\mathrm{}})\right)w_{k,\mathrm{}},\hfill \\ \hfill K^{}w_{k,\mathrm{}}& =\left((\alpha ^k+\alpha ^k)i(\beta ^{\mathrm{}}+\beta ^{\mathrm{}})\right)b_{k,\mathrm{}}.\hfill \end{array}$$

Thus, $`w_{k,\mathrm{}}`$ and $`b_{k,\mathrm{}}`$ are singular vectors and $`\sigma _{k,\mathrm{}}`$ are singular values.

Corollary 5.2:  Let $`𝒢`$ be a $`m\times n`$ rectangular grid. Let

$$\alpha =\mathrm{exp}\left(\frac{\pi i}{m+1}\right),\beta =\mathrm{exp}\left(\frac{\pi i}{n+1}\right),N=\frac{mn}{2}.$$

Then

$$\underset{(k,\mathrm{})X_{m,n}^{}}{}(t(\alpha ^k+\alpha ^k)^2(\beta ^{\mathrm{}}+\beta ^{\mathrm{}})^2)=\underset{j=0\mathrm{}N}{}t^{Nj}(1)^j\underset{||=2j}{}\delta (𝐩_{𝒢,},),$$

where $``$ ranges over all balanced subgraphs with $`2j`$ vertices.

Proof:  This follows directly from Theorem 3.3 and Theorem 4.1.

These results show that the characteristic polynomial of $`KK^{}`$ usually factors a lot if $`𝒢`$ is a rectangle. If $`\zeta `$ is a root of unity whose order $`M`$ is the least common multiple of $`2(m+1)`$ and $`2(n+1)`$ then all the roots $`\sigma _{k,\mathrm{}}^2`$ of this polynomial are in $`RZ[\zeta ]`$, a ring of degree $`\varphi (M)/2`$ over $`Z`$. In particular, for square grids of order $`n`$, irreducible factors of the characteristic polynomial of $`KK^{}`$ have degree at most $`n`$. Here are a few sample examples; we give the polynomial $`det(tIKK^{})=t^n+k_1t^{n1}+\mathrm{}+k_{n1}t+k_n`$ (whose roots are the squares of singular values) factored in $`Z`$.

$$1$$
$`[1,1]`$
$$t1$$
$`[2,1]`$
$$\left(t2\right)^2$$
$`[2,2]`$
$$t2$$
$`[3,1]`$
$$\left(t1\right)\left(t3\right)^2$$
$`[3,2]`$
$$\left(t2\right)^2\left(t4\right)^2$$
$`[3,3]`$
$$t^23t+1$$
$`[4,1]`$
$$\left(t^25t+5\right)^2$$
$`[4,2]`$
$$\left(t^23t+1\right)\left(t^27t+11\right)^2$$
$`[4,3]`$
$$\left(t^26t+4\right)^2\left(t3\right)^4$$
$`[4,4]`$
$$\left(1+t\right)\left(t3\right)$$
$`[5,1]`$
$$\left(t1\right)\left(t2\right)^2\left(t4\right)^2$$
$`[5,2]`$
$$\left(t1\right)\left(t2\right)\left(t5\right)^2\left(t3\right)^3$$
$`[5,3]`$
$$\left(t^23t+1\right)\left(t^25t+5\right)^2\left(t^29t+19\right)^2$$
$`[5,4]`$
$$\left(t1\right)^2\left(t2\right)^2\left(t3\right)^2\left(t6\right)^2\left(t4\right)^4$$
$`[5,5]`$
$$t^35t^2+6t1$$
$`[6,1]`$
$$\left(t^38t^2+19t13\right)^2$$
$`[6,2]`$
$$\left(t^35t^2+6t1\right)\left(t^311t^2+38t41\right)^2$$
$`[6,3]`$
$$\left(t^619t^5+142t^4529t^3+1017t^2922t+281\right)^2$$
$`[6,4]`$
$$\left(t^35t^2+6t1\right)\left(t^38t^2+19t13\right)^2\left(t^314t^2+63t91\right)^2$$
$`[6,5]`$
$$\left(t^310t^2+24t8\right)^2\left(t^310t^2+31t29\right)^4$$
$`[6,6]`$
$$\left(t2\right)^2\left(t^24t+2\right)^2\left(t^28t+8\right)^2\left(t^28t+14\right)^4\left(t4\right)^6$$
$`[7,7]`$
$$(t2)^2(t^312t^2+36t8)^2(t^39t^2+24t17)^4(t^312t^2+45t53)^4$$
$`[8,8]`$
6. Examples

Although Aztec diamonds have so many interesting properties (see \[EKLP\] and \[P\]), the characteristic polynomial of $`KK^{}`$ does not factor very much:

$$(t^38t^2+17t8)^2$$
2-Aztec diamond
$$(t^410t^3+28t^224t+4)^2(t4)^4$$
3-Aztec diamond
$$(t^{10}32t^9+441t^83424t^7+16432t^650240t^5+97041t^4112896t^3+70921t^218784t+1024)^2$$
4-Aztec diamond
$$\genfrac{}{}{0pt}{}{(t^{11}34t^{10}+496t^94064t^8+20562t^766524t^6+137728t^5}{177120t^4+131825t^349066t^2+6576t128)^2(t4)^8}$$
5-Aztec diamond
The fact that these polynomials are always squares follows from symmetry. The factor $`(t4)^{4k}`$ seems to appear in the $`2k+1`$-Aztec diamond, a fact for which we have no explanation.

Finally, here is a small list of “irregular” examples. A possible real Kasteleyn matrix is indicated by the dashed lines (the $`1`$’s).

$$\genfrac{}{}{0pt}{}{t^37t^2+15t9=(t1)(t3)^2}{1.732051,1.732051,1}$$

$$\genfrac{}{}{0pt}{}{t^36t^2+9t4=(t4)(t1)^2}{2,1,1}$$

$$\genfrac{}{}{0pt}{}{t^410t^3+35t^250t+25=(t^25t+5)^2}{1.902113,1.902113,1.175571,1.175571}$$

$$\genfrac{}{}{0pt}{}{t^48t^3+20t^216t+4=(t^24t+2)^2}{1.847759,1.847759,0.765367,0.765367}$$

$$\genfrac{}{}{0pt}{}{{\displaystyle \genfrac{}{}{0pt}{}{t^513t^4+63t^3140t^2+140t49=}{(t^26t+7)(t^37t^2+14t7)}}}{2.101003,1.949856,1.563663,1.259280,0.867768}$$

$$\genfrac{}{}{0pt}{}{{\displaystyle \genfrac{}{}{0pt}{}{t^513t^4+62t^3132t^2+121t36=}{(t4)(t^49t^3+26t^228t+9)}}}{2.126757,2,1.576415,1.197126,0.747468}$$

References:

\[EKLP\] N. Elkies, G. Kuperberg, M. Larsen and J. Propp, Alternating-sign matrices and domino tilings, Journal of Algebraic Combinatorics, 1, 111-132 and 219–234 (1992).

\[K\] P. W. Kasteleyn, The statistics of dimers on a lattice I. The number of dimer arrangements on a quadratic lattice, Phisica 27, 1209-1225 (1961).

\[LL\] E. H. Lieb and M. Loss, Fluxes, Laplacians and Kasteleyn’s theorem, Duke Math. Jour., 71, 337-363 (1993).

\[P\] J. Propp, Twenty open problems in enumeration of matchings, preprint (1996).

\[ST\] N. C. Saldanha and C. Tomei, An overview of domino and lozenge tilings, Resenhas IME-USP, 2(2), 239–252, (1995).

Nicolau C. Saldanha
Departamento de Matemática, PUC-Rio
Rua Marquês de São Vicente 225
Gávea, Rio de Janeiro, RJ 22453-900, BRAZIL
UMPA, ENS-Lyon, 46 Allee d’Italie
69364 Lyon cedex 07, FRANCE
nicolau@mat.puc-rio.br
http://www.mat.puc-rio.br/$``$nicolau/


text:
# Mutual Information of Population Codes and Distance Measures in Probability Space

## Abstract

We studied the mutual information between a stimulus and a system consisting of stochastic, statistically independent elements that respond to a stimulus. Using statistical mechanical methods the properties of the Mutual Information (MI) in the limit of a large system size, N, are calculated. For continuous valued stimuli, the MI increases logarithmically with N and is related to the log of the Fisher Information of the system. For discrete stimuli the MI saturates exponentially with N. We find that the exponent of saturation of the MI is the Chernoff Distance between response probabilities that are induced by different stimuli.

Population coding serves as a common paradigm of information processing in the brain. Its starting point is the fact that very often the response of a single neuron to a stimulus is noisy and is only weakly tuned to changes in the stimulus value. Hence, the information carried by a typical neuron is rather low. The brain may overcome this limitation by distributing the information across a large number of neurons which together carry accurate information about the stimulus. This paradigm has inspired numerous studies of the statistical efficiency of population codes -. Of particular interest is the way in which the amount of information about the stimulus depends on both the size of the responding population and the response properties of the individual neurons.

Many theoretical studies have employed the well-known concept of Fisher Information (FI). The FI is related to the derivative of the population response probability with respect to the stimulus. For a statistically independent population, the FI is an extensive quantity and is relatively easy to calculate. However, it is restricted to the case of a continuously varying stimulus. An alternative measure, applicable for arbitrary stimulus spaces, is provided by Shannon Mutual Information (MI). Unfortunately, except for special cases, exact calculation of the MI for a large population is difficult even for independent populations . This is because MI is bounded from above by the stimulus entropy, and thus is not an extensive quantity. Recently, a relationship between the MI of a continuous stimulus and the FI in a large population has been derived. However, little theoretical progress has been made on the properties of the MI in large systems with discrete stimuli.

In this paper we introduce statistical mechanical methods to calculate analytically the behavior of the MI as the system size, $`N`$, of the population grows. For continuous valued stimuli our theory yields a logarithmic dependence of the MI on the FI, in agreement with previous results . In the case of discrete stimuli, the MI saturates exponentially fast with $`N`$. We show that the exponential saturation rate is dominated by a contribution from the two stimulus values that induce the closest response-probabilities in the population. The contribution of this pair of stimuli to the saturation rate equals the Chernoff Distance between the corresponding response-probabilities, of Large Deviation Theory  .

We consider a population of $`N`$ stochastic units, which we call neurons, that respond simultaneously to a presentation of a stimulus. We denote their responses by a vector $`𝐫=\{r_1,r_2,\mathrm{},r_N\}`$ where $`r_i`$ represents the stochastic response of the i-th neuron to a stimulus. The stimulus states are denoted by the scalar variable $`\theta `$, which can be either discrete or continuous, with a prior probability (or density) $`p(\theta )`$. We denote the probability (or the density) of $`𝐫`$ given a stimulus $`\theta `$ by $`P(𝐫|\theta ).`$ In this paper we will focus on the case of statistically independent neurons (given a stimulus $`\theta `$), namely,

$$P(𝐫|\theta )=\underset{i=1}{\overset{N}{}}P_i(r_i|\theta ).$$
(1)
An important issue is how to quantify the efficiency of the coding of $`\theta `$ in the population responses. For probabilities that are differentiable functions of a continuously varying stimulus, a well-known measure of the efficiency of the population code is the FI. In our case it is

$`J(\theta )`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{^2}{\theta ^2}}\mathrm{ln}P(r_i|\theta )_{r|\theta }`$ (2)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle \frac{}{\theta }}\mathrm{ln}P(r_i|\theta )\right\}^2_{r|\theta }`$ (3)
where $`\mathrm{}_{r|\theta }`$ denotes an average with respect to $`P(r_i|\theta )`$, which is clearly an extensive quantity. Here we study an alternative measure of efficiency, the MI of the system, which is the average amount of information on $`\theta `$ that is added by observing the response, $`𝐫`$. It is useful to define a local MI, $`I(\theta )`$,

$$I(\theta )\mathrm{log}_2\left\{\frac{P(𝐫|\theta )}{P(𝐫)}\right\}_{𝐫|\theta }$$
(4)
where $`P(𝐫)=_\theta P(𝐫|\theta )p(\theta )`$, and $`\mathrm{}_{𝐫|\theta }`$ denotes an average w.r.t. $`P(𝐫|\theta ).`$ The full MI, $`I`$, is

$$I=I(\theta )_\theta $$
(5)
where $`\mathrm{}_\theta `$ is an average over the stimulus distribution. A central quantity is the log-likelihood variable,

$$S(𝐫,\varphi ,\theta )=\mathrm{ln}\left[\frac{P(𝐫|\varphi )}{P(𝐫|\theta )}\right]=\underset{i=1}{\overset{N}{}}S_i(r_i,\varphi ,\theta )$$
(6)
where $`S_i(r_i,\varphi ,\theta )=\mathrm{ln}P(r_i|\varphi )\mathrm{ln}P(r_i|\theta )`$. Equation (4) can be written as

$$I(\theta )=\underset{\varphi \theta }{}dS_\varphi P_\theta \{S_\varphi \}\mathrm{log}_2\mathrm{exp}S_\varphi ^{}_\varphi ^{}$$
(7)
where $`P_\theta \{S_\varphi \}`$ is the distribution of $`S_\varphi `$ $`S(𝐫,\varphi ,\theta )`$, calculated with respect to $`P(𝐫|\theta ).`$ For large $`N`$ this distribution is centered around the mean value of $`S_\varphi `$

$$\underset{i=1}{\overset{N}{}}S_i(r_i,\varphi ,\theta )_{r|\theta }=D_{KL}(\varphi ||\theta )$$
(8)
where $`D_{KL}(\varphi ||\theta )>0`$ is the relative entropy of $`P(𝐫|\varphi )`$ and $`P(𝐫|\theta ),`$ also known as the Kullback-Leibler (KL) distance between the two distributions. The correlation matrix of the fluctuations $`\delta S_\varphi =`$ $`S_\varphi +D_{KL}(\varphi ||\theta )`$ is also of order $`N`$, and is given by

$$C(\varphi ,\psi |\theta )=\underset{i=1}{\overset{N}{}}\delta S_i(r_i,\varphi ,\theta )\delta S_i(r_i,\psi ,\theta )_{r|\theta }.$$
(9)
We first discuss the MI of a large system with a continuous varying stimuli. Here, all averages over stimulus space stand for integrals with a density $`p(\theta )`$, and we assume that the probabilities vary smoothly with $`\theta .`$ In this case, the dominant contribution to Eq. (7) comes from values of $`S_\varphi `$ near their mean value and $`\varphi `$ which is near $`\theta `$, i.e., $`|\varphi \theta |1/\sqrt{N}.`$ This is because $`S_\theta =S(𝐫,\theta ,\theta )=0`$ while for $`\varphi `$ that is far from $`\theta `$, $`S_\varphi D_{KL}(\varphi ||\theta )`$, which is large and negative. For small magnitudes of $`\delta \varphi `$ =$`\varphi \theta `$, and $`\delta \psi `$ =$`\psi \theta `$ we can write

$$D_{KL}(\varphi ||\theta )\frac{1}{2}J(\theta )\delta \varphi ^2,C(\varphi ,\psi |\theta )J(\theta )\delta \varphi \delta \psi $$
(10)
where $`J(\theta )`$ is the FI, Eq. (2). The low rank form of the correlation matrix in Eq. (10) implies that the fluctuations $`\delta S_\varphi `$ can be described as $`z\sqrt{J(\theta )}\delta \varphi `$, where by central limit theorem $`z`$ is a Gaussian variable with zero mean and unit variance. Substituting this in Eq. (7) yields,

$`I(\theta ){\displaystyle }Dz\times `$ (11)
$`\mathrm{log}_2\left\{{\displaystyle 𝑑\psi p(\psi +\theta )\mathrm{exp}\left(\frac{J(\theta )\psi ^2}{2}+z\sqrt{J(\theta )}\psi \right)}\right\}`$ (12)
where $`Dz=dz\mathrm{exp}(z^2/2)/\sqrt{2\pi }.`$ Evaluating the integrals in the limit of a large $`J`$, and substituting in Eq. (5), yields

$$IH_\theta \frac{1}{2}\mathrm{log}_2\left\{\frac{2\pi e}{J(\theta )}\right\}_\theta $$
(13)
where $`H_\theta =𝑑\theta p(\theta )\mathrm{log}_2p(\theta )`$, in agreement with previous results.

We now turn to the more difficult case of discrete valued stimulus which takes the values $`\left\{\theta _l\right\}_{l=1}^M`$ and $`M`$ remains finite as $`N\mathrm{}`$. The term $`\varphi =\theta `$ in the average $`\mathrm{exp}S_\varphi _\varphi `$ in Eq. (7) is $`p(\theta )`$, yielding a total contribution to $`I`$ which is the stimulus entropy $`H_\theta =_\theta p(\theta )\mathrm{log}_2p(\theta )`$. Naively, we would therefore expect that the main contribution to $`\mathrm{log}_2p(\theta )I(\theta )`$ comes from the typical values of $`S_\varphi `$, namely, $`D_{KL}(\varphi ||\theta )`$, for a state $`\varphi `$ which is closest to $`\theta .`$ Such a contribution would be proportional to $`\mathrm{exp}(D_{KL}(\varphi ||\theta ))`$. However, we find that in fact, the dominant contribution comes from rare values of $`𝐫`$ such that $`S_\varphi `$ $`0`$ for one of the states $`\varphi `$ $`(\varphi \theta )`$. This is because, as we will show, although this regime has an exponentially small probability, its contribution to $`I`$ is exponentially larger than that of the typical value of $`S_\varphi `$ making it the dominant correction to $`IH_\theta `$.

To evaluate Eq. (7) in the discrete case we use an integral representation of the distribution of $`S_\varphi `$ $`S(𝐫,\varphi ,\theta )`$,

$`P_\theta \{S_\varphi \}={\displaystyle \underset{\varphi \theta }{}\frac{dY_\varphi }{2\pi }\mathrm{exp}\{F_\theta (Y_\varphi ,S_\varphi )\}}`$ (14)
$`F_\theta ={\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{ln}\mathrm{exp}i\mathrm{\Sigma }_\varphi Y_\varphi S_i(r_i,\varphi ,\theta )_{r|\theta }+i{\displaystyle \underset{\varphi }{}}Y_\varphi S_\varphi .`$ (15)
Note that there is no integration over variables with $`\varphi =\theta `$, since $`S_\theta `$ $`S(𝐫,\theta ,\theta )`$ $`=0`$. The large $`N`$ limit of Eq. (7) and (14) is evaluated by the saddle point method. Solving the saddle point equations for $`\{S_\varphi \}`$, and $`\{Y_\varphi \}`$ we find that there are $`M1`$ saddle points, each of which is characterized by having one $`S_\varphi `$ of $`O(1)`$ while the remaining $`S_\varphi ^{}`$ are negative and of $`O(N)`$. The auxiliary variable $`Y_\varphi =i\alpha `$ where $`\alpha `$ is a real number of order $`1`$ while the remaining $`Y_\varphi ^{}`$ are zero. At this saddle point, the value of $`F_\theta `$, Eq. (15) is

$$D_\alpha (\varphi ||\theta )\underset{i=1}{\overset{N}{}}\mathrm{ln}\mathrm{exp}\alpha S_i(r_i,\varphi ,\theta )_{r|\theta }$$
(16)
where the order parameter $`\alpha `$ is evaluated by maximizing Eq. (16),

$$\frac{}{\alpha }D_\alpha (\varphi ||\theta )=0.$$
(17)
We will denote the maximum value of $`D_\alpha `$ by $`D_C(\varphi ,\theta )`$. Out of these $`M1`$ saddle points the dominant one is that with the smallest $`D_C(\varphi ,\theta )`$, yielding

$$\mathrm{ln}\left[\mathrm{log}_2p(\theta )I(\theta )\right]\underset{\varphi \theta }{\mathrm{min}}D_C(\varphi ,\theta ).$$
(18)
Finally, the full MI is given by

$$\mathrm{ln}\left[H_\theta I\right]=\underset{\left(\varphi \psi \right)}{\mathrm{min}}D_C(\varphi ,\psi )\frac{1}{2}\mathrm{ln}N+A$$
(19)
where the $`\mathrm{ln}N`$ and the constant $`A`$ correction terms come from the Gaussian integration around the saddle point of Eqs. (7) and (14). The calculations of the constant $`A`$ will be presented elsewhere . Finally, examining Eq. (14) at the saddle-point, it can be seen that $`\mathrm{exp}(D_C(\varphi ,\theta ))`$ equals the probability that $`S_\varphi /N0`$ for a large $`N`$. This implies, as we have mentioned above, that the rate of saturation of the MI is dominated by the probability of the rare event that one of the log-likelihood ratios is close to zero.

The quantity $`D_C(\varphi ,\psi )`$ derived in the above theory is identical to the Chernoff Distance (or Chernoff Information) between the two distributions $`P(𝐫|\varphi )`$ and $`P(𝐫|\psi )`$. More generally, the Chernoff Distance between two arbitrary distributions $`P(x)`$, and $`Q(x)`$, is defined  as

$$D_C(Q,P)=\underset{\alpha }{\mathrm{max}}D_\alpha (Q||P)$$
(20)
where

$$D_\alpha (Q||P)\mathrm{ln}\underset{x}{}Q(x)^\alpha P(x)^{1\alpha }$$
(21)
$`D_\alpha (Q||P)`$ are proportional to the family of Renyi distances.. $`D_\alpha (Q||P)`$ vanishes at $`\alpha =0`$ and $`\alpha =1`$. It is positive (if $`QP`$) for $`0<\alpha <1`$, with a maximum at $`\alpha ^{}`$, $`0<\alpha ^{}<1`$, and is negative outside this regime. It is related to the KL distance through its slope at $`\alpha =0`$, i.e., $`D_\alpha (Q||P)/\alpha |_{\alpha =0}=D_{KL}(Q||P)`$. The Chernoff Distance chooses the value of $`\alpha `$ which maximizes $`D_\alpha (Q||P)`$. This value is not constant but depends on the pair of distributions $`Q`$ and $`P`$. Note, that in the case of a family of distributions parameterized by a continuous parameter $`\theta `$, the Chernoff Distance is related to the FI through $`4D_C(\theta ,\theta +\delta \theta )D_{KL}(\theta ,\theta +\delta \theta )J(\theta )(\delta \theta )^2/2`$.

Although $`D_C(Q,P)`$ is not symmetric with respect to $`Q`$ and $`P`$ for general $`\alpha `$, it obeys the symmetry $`D_\alpha (Q||P)=D_{1\alpha }(P||Q),`$ which implies that $`D_\alpha ^{}(Q||P)=`$ $`D_\alpha ^{}(P||Q)`$. Thus, the Chernoff Distance of $`Q`$ and $`P`$ is a symmetric function of the two distributions. It is smaller than $`D_{KL}`$, and in addition, it is less sensitive than the KL distance to outlier states. In particular, a single state which has a nonzero probability $`P`$ but zero probability $`Q`$, causes $`D_{KL}(Q||P)`$ to diverge whereas $`D_C(Q,P)`$ remains finite. In fact, it diverges only when the two distributions have zero overlap, i.e., the intersection of their support is empty.

Equation (19) implies that the Chernoff distance controls the rate of saturation of the MI. In order to test our theory we have studied the case of a population of statistically independent, binary neurons, where $`r_i=\{0,1\},`$ and a stimulus which can take three values: $`\theta _1,\theta _2`$ and $`\theta _3`$. Each neuron has a preferred stimulus, denoted by $`\theta ^i.`$ The mean value of $`r_i`$ (i.e., the probability that $`r_i=1`$) is

$$f_i(\theta )=z_i(\theta )+T\delta _{\theta ,\theta ^i}$$
(22)
For each neuron, the $`\theta _i`$ is chosen at random from $`\{\theta _1,\theta _2,\theta _3\},`$ with equal probability, and $`z_i(\theta )`$ are chosen at random uniformly from $`[a,b]`$, with $`Tba`$. Thus, the parameter $`T`$ measures the selectivity of the responses to the stimulus values. In this case, there is no statistical difference between the response probability of the population as a whole to the three stimuli, and changing the stimulus corresponds to permuting the mean responses of the individual neurons in the population. This is common in biological situations where different stimuli (such as different angles, spatial positions or abstract objects) elicit activity profiles that are similar in shape but shifted in position across the network.

Figure 1 shows a nice agreement between the results of exact numerical calculation of the MI of this model up to size $`N=25,`$ and the asymptotic result of Eqs. (16)-(19), evaluated for the above distributions. Note that because of the symmetry between the different stimuli in this example, $`D_C(\varphi ,\psi ),`$ is the same for all pairs of stimuli, $`\theta \theta ^{}`$, so that they all contribute equally to Eq. (19).

The importance of $`D_C(\varphi ,\psi )`$ in characterizing the efficiency of the population code is manifest not only in the MI of the system but also in the accuracy of the discrimination between stimuli on the basis of observation of the population responses. Plausible discriminatory often base their discrimination between a pair of stimuli on the log-likelihood ratio of the corresponding distributions. In particular, the Maximum-Likelihood (ML) discriminator makes a deterministic decision between a pair of stimuli $`\theta `$ and $`\varphi `$ upon observing $`𝐫`$, according to whether the log-likelihood, $`S(𝐫,\varphi ,\theta )`$, is larger or smaller than zero. As outlined below, our theory yields that the probability of confusion of an ML discrimination between a pair of stimuli, $`\theta `$ and $`\varphi ,`$ is dominated by the probability that $`S(𝐫,\varphi ,\theta )`$ is close to zero, and thus is determined by $`D_C(\varphi ,\theta )`$.

To evaluate the discrimination error, we write the average probability of confusion, $`P_C`$ of an ML discrimination as $`P_C=P_C(\theta )_\theta `$ where,

$$P_C(\theta )=1\underset{\varphi \theta }{}dS_\varphi \mathrm{\Theta }(S_\varphi )P_\theta \{S_\varphi \}\underset{\varphi \theta }{}$$
(23)
where $`\mathrm{\Theta }(x)`$ is the Heaviside step function, and $`P_\theta \{S_\varphi \}`$ is given by Eq. (14). Using the saddle point methods as described above, we find that in a large system the dominant contribution to the integrals over $`S_\varphi `$ comes from the edge points where one $`S_\varphi `$ is close to zero. Evaluating the saddle point equations for $`Y_\varphi `$ under this condition, yields

$$\mathrm{ln}P_C=\underset{\left(\varphi \psi \right)}{\mathrm{min}}D_C(\varphi ,\psi )\frac{1}{2}\mathrm{ln}N+A^{}$$
(24)
An example is shown in Fig. 2, where we have computed the confusion error of an ML discriminator between two stimuli in the case of a population of neurons responding to three stimuli as described above, except that in this case, the neurons are Poissonian with means $`f_i(\theta )`$, of the form Eq. (22). The numerical results obtained by simulating the ML discriminator are in very good agreement with the prediction of Eq. (24).

In conclusion, we have shown the relation between the MI of a large population coding for a stimulus and the distance between the response-probabilities that are induced by the different stimulus values. In the case of a continuous stimulus, the MI increases logarithmically with $`N`$, for large $`N`$, and is related to the FI which measures the vanishing rate of the distance for infinitesimally small stimulus differences. In the case of discrete stimulus, the MI saturates exponentially at a rate which is given by the Chernoff Distance between the closest pair of population response probabilities. In addition, we have shown that $`D_C`$ determines also the probability of discrimination error. This extends the classical Large Deviation Theory results to cases where the elements of the population are not identical. Our finding that the Chernoff Distance controls both the MI and the error probability for a large population that code discrete stimuli, suggests that $`D_C`$ is a useful measure of the quality of neuronal population codes. We hope that these results will provide tractable tools to study the nature of population codes in the brain using experimental data on neuronal representations of sensory, motor and cognitive events.

This work was partially supported by grants of the Israeli Science Foundation and the Israel-USA Binational Science Foundation. We thank Samengo and Treves for discussions on Ref..


text:
# Total Photonic and Hadronic Cross-sections

## I Introduction

This talk will be a short review of the status of the calculation of total cross-sections using a QCD driven mini-jet modelemmus , with particular emphasis on recent measurements and results of theoretical calculations of $`\gamma \gamma `$ cross-sections.

One of the aim of QCD is to calculate cross-sections for hadronic processes. We find that the level of available experimental information on total hadronic cross-sections has now reached a stage so as to allow, for the first time, definite progress in the calculation of the one quantity which has so far escaped a complete quantitative understanding in a QCD framework, namely the total hadronic cross-section. We have a complete set of processes, $`pp`$, $`p\overline{p}`$, $`\gamma p`$ and $`\gamma \gamma `$ measured in a common energy range, $`\sqrt{s}=1÷100GeV`$, with the purely hadronic processes measured up to $`\sqrt{s}3\times 10^5GeV`$ martin . The latter allows for very good parametrizations in a large energy range, the other two cross-sections to test the QCD content of hadrons versus the one in the photons. In Figure 1 we show a compilation of all presently available photon and proton total cross-sections, scaling them so as to be able to compare one to each other.

For $`\gamma p`$ processes, the scale factor is the product of a Quark Parton Model factor $`3/2`$ multiplied by a Vector Meson Dominance type factor $`1/240`$emmus , for $`\gamma \gamma `$ we just square this factor.

The plan of this paper is as follows :

1. In Sect. II we describe our recent results for $`\gamma \gamma `$ and extrapolate them to $`\gamma p`$
2. In sect. III we discuss a possible ansätz on an extension of the minijet model to $`\gamma ^{}\gamma ^{}`$
3. In Sect. IV we present recent results for $`pp`$ and $`p\overline{p}`$, using an impact parameter distribution for partons in the hadrons obtained from the Bloch-Nordsieck summation technique, and discuss its possible extensions to real and virtual photons.

## II $`\gamma \gamma `$ and $`\gamma proton`$

While for quite some time the photon-photon data at LEP exhibited a discrepancy between different collaborations (although within the experimental errors), we now have two sets of data points which are in excellent agreement with each other. A theorist can then start his/her work. We show in Figure 2 the description of L3L3  and OPAL dataOPAL  using the Eikonal Minijet Model (EMM).

The theoretical context in which this curve was obtained is discussed in emmus ; zeitus  : the EMM uses the eikonal approximation eikminijets  to calculate the total cross-section, i.e.

$$\sigma _{tot}=2P_{had}^{ab}d^2\stackrel{}{b}[1e^{i\chi (b,s)}]$$
(1)
and approximates the eikonal function neglecting $`\mathrm{}e\chi (b,s)`$ and putting

$$2\mathrm{}m\chi (b,s)=n(b,s)=A(b)[\sigma _{soft}(s)+\frac{\sigma _{jet}(s,p_{tmin})}{P_{had}^{ab}}]$$
(2)
where the QCD calculable jet cross-section is the quantity which drives the rise therise  in all total cross-sections. This function is defined as

$$\sigma _{jet}=d^2\stackrel{}{p}_t\frac{d\sigma ^{QCD}}{d^2\stackrel{}{p}_t}=𝑑x_1𝑑x_2f^{i/a}(x_1)f^{j/b}(x_2)d^2\stackrel{}{p}_t\frac{d\widehat{\sigma }^{ij,kl}}{d^2\stackrel{}{p}_t}$$
(3)
where $`\widehat{\sigma }`$ is the parton-parton cross-section for the subprocess $`ijkl`$, the sum runs over i,j,k,l,=parton type and the integration covers a region from a minimum $`p_t`$ to the entire phase space. The quantity $`\sigma _{jet}`$ is higly dependent upon the regulator $`p_{tmin}`$ and one of the aims of a complete QCD calculation is to eliminate this dependence. Presently this is not yet possible, but one can nonetheless expect $`p_{tmin}`$ to be the smallest momentum exchanged between partons such that perturbative QCD can be applied, namely not less than 1 GeV and probably not more than 2 GeV. The parameter $`P_{had}^{ab}P_{had}^aP_{had}^b`$ is in principle energy dependent and $`P_{had}^a`$ can be interpreted as the probability for particle $`a`$ to behave like a hadron, i.e. $`P_{had}^{hadron}1`$ and, typically, $`P_{had}^\gamma 𝒪(\alpha )`$ ladinsky ; halzen . Phenomenologically, in order to obtain a description of the measured total cross-sections, there are other input quantities which need to be fixed, namely $`\sigma _{soft}`$ and the b-distribution function $`A(b)`$. One way to proceed has been, previously, to determine all the parameters from the process $`\gamma p`$. However, the EMM model is not working very well for the proton case, as discussed in ff2 , probably because of uncertainties in the hadronic transverse momentum distributions, and we have opted, in this note, for a different approach, namely we obtain the parameters from $`\gamma p`$, extrapolate them to $`\gamma \gamma `$ varying them within at most 10% to fit the $`\gamma \gamma `$ cross-section, and then revert back and see what the best description of $`\gamma \gamma `$ will produce, when applied to photoproduction. The result of this modus operandi is presented here. For Figure 2, we have chosen to describe the partonic matter distribution inside the hadrons through a function inspired by the pion electromagnetic form factor, namely

$$A(b)=\frac{1}{(2\pi )^2}d^2\stackrel{}{q}e^{i\stackrel{}{q}\stackrel{}{b}}(\frac{k_0^2}{q^2+k_0^2})^2$$
(4)
The scale $`k_0`$ has been let to vary, according to an intrinsic transverse momentum ansätz emmus . Thus, Figure 2 corresponds to $`k_0=0.4GeV`$, $`P_{had}^\gamma =1/240,\sigma _{soft}=(21+\frac{42}{s})mb`$, $`p_{tmin}=1.5GeV`$ and GRS GRS  type densities for the photon. Next we extrapolate this curve to $`\gamma p`$ processes, putting $`\sigma _{soft}^{\gamma p}=3/2\sigma _{soft}`$, the proton form factor with dipole type expression instead of one of the photon type monopole expressions, GRV type densities for the proton GRV  and GRS for the photon. The result is the upper curve shown in Figure 3 and compared with old and recent data HERAZ ; HERAH1 ; DIS ; camil . For completeness and as a reference to our previous work, we also show a band, with the lowest curve corresponding to GRV densities for both proton and photon, $`p_{tmin}=2GeV`$, $`k_0=0.66GeV`$ for the photon, same $`P_{had}^\gamma `$ as the upper curve, and $`\sigma _{soft}^{\gamma p}=(31+\frac{10}{\sqrt{s}}+\frac{38}{s})mb`$.

## III The case of the virtual photon

Hadronic interactions of the virtual photon in $`e^+e^{}`$ and $`\gamma \gamma `$ processes have been a subject of interest for some time now virt . The appearance of experimental data on jet production with virtual photons at HERA virzeus  and total $`\gamma ^{}\gamma ^{}`$ cross-sections at LEP l3virt , has given added impetus to develop a model which will describe the total cross-sections for virtual photons, especially the ’resolved’ part torbjorn .

We want to develop a model to understand these cross-sections in the context of EMM. We propose that the virtual photon description is the same as for the case of real photons except that the intrinsic transverse momentum ansätz is complemented with a factor inspired by Extended Vector Meson Dominancebernd . At high $`Q_\gamma ^2`$, a factor like

$$A((b,Q_1^2,Q_2^2)=\frac{m_\rho ^2}{m_\rho ^2+Q_1^2}\frac{m_\rho ^2}{m_\rho ^2+Q_2^2}A(b)$$
(5)
will suppress the part of the cross-section which has a hadronic content. In Eq.(5), $`P_{had}^\gamma (Q_\gamma ^2=0)`$ is the usual factor, $`1/240`$, obtained from Vector meson Dominance and discussed in many paperszeitus ; ladinsky ; halzen . In Figure 4 we show our predictions for the hadronic content of the $`\gamma ^{}\gamma ^{}`$ cross-section, excluding for the time being the direct and single resolved contribution. Similar results can also be obtained using factorization martinstar  in the context of the Aspenaspen  model.

## IV The Bloch Nordsieck Model for the impact parameter distribution

It is evident from the previous discussion, and we have pointed this out in many papers, that it is not possible to have a QCD description of the total cross-section without understanding the transverse momentum distribution of partons in the colliding hadrons, or, in the eikonal language, without understanding the impact parameter distribution. We have attempted such description, and will show in the following the main highlights. Some very recent results in the case of proton proton and proton antiproton are shown in Figure 5,

where the two curves have been obtained using the eikonal formula of Eq.(1), the jet cross-sections with GRV densities and $`p_{tmin}=1.15GeV`$, $`\sigma _{soft}^{pp}=48mb`$, $`\sigma _{soft}^{p\overline{p}}=48(1+\frac{2}{\sqrt{s}})mb`$. Unlike the curves obtainable in a straightforward application of the EMM with b-distribution determined by the proton form factor, this figure uses an s-dependent b-distribution given by

$$A(b,s)=\frac{d^2\stackrel{}{q}e^{i\stackrel{}{b}\stackrel{}{q}}\mathrm{\Pi }(q,s)}{(2\pi )^2\mathrm{\Pi }(0,s)}=\frac{e^{h(b,s)}}{d^2\stackrel{}{b}e^{h(b,s)}}$$
(6)
where $`\mathrm{\Pi }(q,s)`$ is the transverse momentum distribution of colliding partons generated by soft gluon radiation from the initial state and, to leading order,

$$h(b,s)=\frac{8}{3\pi }_0^{q_{max}}\alpha _s(k_t)\frac{dk_t}{k_t}[1J_0(bk_t)]\mathrm{log}\frac{q_{max}+\sqrt{q_{max}^2k_t^2}}{q_{max}\sqrt{q_{max}^2k_t^2}}$$
(7)
The functions $`h(b,s)`$ and $`A(b,s)`$ are energy dependent inasmuch as the kinematics of a given subprocess determine how much energy is available to single soft gluon emission. The energy dependence appears through the quantity $`q_{max}`$ which is a function of the c.m. energy $`\sqrt{s}`$ and $`p_{tmin}`$. In addition $`h(b,s)`$ depends upon the single soft gluon distribution function $`\alpha _s(k_t)dk_t/k_t`$. Since $`\mathrm{\Pi }(q,s)`$ is summed over all possible gluon distributions, the infrared divergence is of course cancelled, there remains, however, the problem of evaluating $`\alpha _s`$ in the infrared limit. We have used a phenomenological ansätz, namely a singular but integrable $`\alpha _s`$ and the result shown in Figure 5 is discussed in ff2  for the rising part of the total proton cross-sections. The main characteristics of this treatment is that, as the minijet cross-section rises with energy, soft gluon emission produces an acollinearity of the partons and reduces the probability of collisions. This affects the cross-sections in two ways : at low energy it produces a very soft decrease in $`\sigma ^{pp}`$ and contributes to the faster decrease in $`\sigma ^{p\overline{p}}`$, at high energy it tames the rise due to $`\sigma ^{jet}`$. It is then possible to have a very small $`p_{tmin}`$ to see the very beginning of the rise around $`10÷20GeV`$, without having too large a cross-section when the energy climbs into the TeV range and beyond.

The above discussion illustrates also the way to study virtual and real photon interactions : the aim is to obtain an energy and momentum dependence in the impact parameter distribution through the kinematics characterizing real or virtual photon process, as it has been done for the proton. At high $`Q_\gamma ^2`$, one can expect $`q_{max}`$ in Eq.(7) to increase with $`Q_\gamma ^2`$, thus producing the same suppression effect that high c.m. energy values have in the softening of the rise due to the jets. The overall effect may be similar to the EVMD factor discussed in Sect. III. Work is in progress to apply the Bloch-Nordsieck formalism to these collisions.

## V conclusion

We have discussed total cross-sections for protons and photons, real and virtual, using QCD calculable mini-jets cross-sections as the physical process which drives the rise of all total cross-sections. We have used the eikonal formalism to unitarize the cross-section and discussed possible ways to reduce the arbitrariness introduced in this formalism by the impact parameter distribution. A QCD model using soft gluon summation has been presented and compared with data for purely proton processes.

## VI acknowledgement

Two of us, A.G. and G.P., acknowledge the support of EEC TMR-CT98-0169.