June 01, 2026
We prove that the bilinear generating function for Wronskians of Hermite polynomials can be expressed as the classical Mehler kernel multiplied by a polynomial, thereby extending the result of Pupasov-Maksimov [1] for exceptional Hermite polynomials. We establish several properties of the polynomials appearing in this extended version of the Mehler formula and present four conjectures about them.
: Wronskian, Hermite polynomials, Mehler formula, heat kernel, Darboux transformation
: Primary 26C05, Secondary 33C45
Let \(H_n(x)\) denote the Hermite polynomials. It is known (see [2]) that for all \(|z|<1\) and \(x,y\in {\mathbb{C}}\) \[E(z,x,y):=\sum\limits_{n\geqslant 0} z^n \frac{H_n(x) H_n(y)}{2^n n!}= \frac{1}{\sqrt{1-z^2}} \exp\bigg( \frac{2xzy-z^2(x^2+y^2)}{(1-z^2)} \bigg).\] The expression on the right-hand side is called the Mehler kernel, and it plays a fundamental role in many applications. For example, it is used to construct the kernel of the fractional Fourier transform [3], the propagator for the harmonic oscillator in quantum mechanics [1], and the transition probability density of the Ornstein-Uhlenbeck process [4].
Before stating our results, we introduce some notation. We denote by \({\mathbb{N}}\) the set of positive integers and by \({\mathbb{Z}}_{\geqslant 0}\) the set of nonnegative integers. The Wronskian of smooth functions \(\{f_j(x)\}_{1\leqslant j \leqslant m}\) is defined as \[\textrm{Wr}[f_1,\dots,f_m]:=\det\,[\partial_x^{i-1} f_j(x)]_{1\leqslant i,j \leqslant m}.\] For an \(m\)-tuple \(\boldsymbol{\textsl{j}}=(j_1,j_2,\dots,j_m) \in {\mathbb{Z}}_{\geqslant 0}^m\) we denote \[H_{\boldsymbol{\textsl{j}}}(x):=\textrm{Wr}[H_{j_1}(x), \dots,H_{j_m}(x)].\] We set \(H_{\boldsymbol{\textsl{j}}}(x) \equiv 1\) if \({\mathbf{j}}=\varnothing\) (the empty tuple). In what follows, \(\boldsymbol{\textsl{k}}=(k_1,k_2,\dots,k_l)\in {\mathbb{N}}^l\) will always denote a strictly increasing sequence of positive integers, and \(l\) will denote the length of the sequence \(\boldsymbol{\textsl{k}}\). For \(n\in {\mathbb{Z}}_{\geqslant 0}\) we denote \(H_{\boldsymbol{\textsl{k}},n}(x):=H_{\boldsymbol{\textsl{j}}}(x)\) where \(\boldsymbol{\textsl{j}}=(k_1,k_2,\dots,k_l,n)\).
Wronskians of Hermite polynomials are important objects that arise in several contexts. They play a fundamental role in the classification of rational monodromy-free Schrödinger operators with rational potentials having quadratic growth at infinity [5]. They appear as rational solutions of the fourth Painlevé equation [6]. Zeros of Wronskians of Hermite polynomials were studied in [7]. It is known [8] that \(H_{\boldsymbol{\textsl{k}}}(x)\) has no real zeros if and only if \(\boldsymbol{\textsl{k}}\) is a Krein-Adler sequence, that is, a strictly increasing sequence \(\boldsymbol{\textsl{k}}=(k_1,k_2,\dots, k_l) \in {\mathbb{N}}^l\) such that \(l\) is even and \(k_{2i}=k_{2i-1}+1\) for all \(1\leqslant i \leqslant l/2\).
We now define our main object of interest: \[\label{def:Ek} E_{\boldsymbol{\textsl{k}}}(z,x,y):=\sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} z^n \frac{H_{\boldsymbol{\textsl{k}},n}(x) H_{\boldsymbol{\textsl{k}},n}(y)}{2^{n+l} n! \prod\limits_{j=1}^l (n-k_j)}.\tag{1}\] As we will show later in Lemma 1, the series in 1 converges for \(|z|<1\), \(x,y\in {\mathbb{C}}\) and defines an analytic function of three variables in that domain. Our main result is the following
Theorem 1. For any strictly increasing sequence \(\boldsymbol{\textsl{k}}\in {\mathbb{N}}^l\) there exist polynomials \(Q^{\boldsymbol{\textsl{k}}}_n(x,y)\), \(0\leqslant n \leqslant k_l+1\) such that \[\label{eqn95main} E_{\boldsymbol{\textsl{k}}}(z,x,y)=E(z,x,y) \sum\limits_{n=0}^{k_l+1} z^n Q^{\boldsymbol{\textsl{k}}}_n(x,y),\qquad{(1)}\] for all \(|z|<1\), \(x,y\in {\mathbb{C}}\).
Theorem 1 was established by A. M. Pupasov-Maksimov [1] in the important case where \(\boldsymbol{\textsl{k}}\) is a Krein-Adler sequence. To explain the significance of this result and the role that Krein-Adler sequences play, let us consider the functions \[\label{def95fkn} f_{\boldsymbol{\textsl{k}},n}(x):=e^{-x^2/2} \frac{H_{\boldsymbol{\textsl{k}},n}(x)}{H_{\boldsymbol{\textsl{k}}}(x)}.\tag{2}\] It is known that for every strictly increasing sequence \(\boldsymbol{\textsl{k}}\in {\mathbb{N}}^l\) these functions are solutions of the second-order linear differential equation \[\label{f95kn95ODE} {\mathcal{L}} f(x):=-f''(x)+U_{\boldsymbol{\textsl{k}}}(x) f(x)=(2n+1-2l) f(x),\tag{3}\] where \(U_{\boldsymbol{\textsl{k}}}(x)\) is a rational function defined by \[U_{\boldsymbol{\textsl{k}}}(x):= x^2- 2\partial_x^2 \ln(H_{\boldsymbol{\textsl{k}}}(x))= x^2 - 2\frac{H_{\boldsymbol{\textsl{k}}}''(x)}{H_{\boldsymbol{\textsl{k}}}(x)}+ 2\frac{H_{\boldsymbol{\textsl{k}}}'(x)^2}{H_{\boldsymbol{\textsl{k}}}(x)^2}.\] This result follows from [9] and it is stated in this exact form in the proof of [9]. The differential operator in 3 and the solutions of equation 3 can be obtained by a sequence of Darboux transformations, see [5] and [9].
As mentioned above, when \(\boldsymbol{\textsl{k}}\) is a Krein-Adler sequence, the polynomial \(H_{\boldsymbol{\textsl{k}}}(x)\) has no real zeros (see [9]) and the polynomials \(H_{\boldsymbol{\textsl{k}},n}(x)\) are called exceptional Hermite polynomials. In this case, the functions \(f_{\boldsymbol{\textsl{k}},n}(x)\) belong to the Schwartz class (the class of smooth rapidly decreasing functions on \({\mathbb{R}}\)) and satisfy the orthogonality condition \[\int_{{\mathbb{R}}} f_{\boldsymbol{\textsl{k}},n}(x) f_{\boldsymbol{\textsl{k}},m}(x) {\textrm d}x= \delta_{n,m} \sqrt{\pi} 2^{n+l} n! \prod_{j=1}^l (n-k_j), \;\;\; n,m \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}.\] Moreover, the functions \(\{f_{\boldsymbol{\textsl{k}},n}(x)\}_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}}\) form a complete orthogonal basis of \(L_2({\mathbb{R}}, {\textrm d}x)\). These results can be found in [9], [10]. Thus, when \(\boldsymbol{\textsl{k}}\) is a Krein-Adler sequence, the functions \(\{f_{\boldsymbol{\textsl{k}},n}(x)\}_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}}\) form a complete eigenbasis for the operator \({\mathcal{L}}\) in \(L_2({\mathbb{R}}, {\textrm d}x)\). This allows one to write down the spectral expansion of the corresponding heat kernel [11], [12] (the integral kernel of the operator \(\exp(-t {\mathcal{L}})\)) in the form \[\label{L95heat95kernel} p(t,x,y)=\sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} e^{-(2n+1-2l) t} \; \frac{f_{\boldsymbol{\textsl{k}},n}(x) f_{\boldsymbol{\textsl{k}},n}(y)}{\lVert f_{\boldsymbol{\textsl{k}},n} \rVert^2}=\frac{1}{\sqrt{\pi}} e^{(2l-1)t-(x^2+y^2)/2} \;\frac{E_{\boldsymbol{\textsl{k}}}(e^{-2t},x,y)}{H_{\boldsymbol{\textsl{k}}}(x) H_{\boldsymbol{\textsl{k}}}(y)}.\tag{4}\] Thus, the significance of Theorem 1 in the Krein-Adler case is that it gives an explicit expression for the heat kernel of the Schrödinger operator \({\mathcal{L}}\).
Remark 1: Pupasov-Maksimov [1] states his results using the probabilists’ version of Hermite polynomials, defined by \({\textrm{He}}_n(x):=2^{-n/2} H_n(x/\sqrt{2})\). Denoting the polynomials appearing in [1] by \(\widetilde{Q}_n^{\boldsymbol{\textsl{k}}}(x,y)\), the relation between these two families of polynomials is \[\widetilde{Q}_n^{\boldsymbol{\textsl{k}}}(x,y)=C \times Q_n^{\boldsymbol{\textsl{k}}}(x/\sqrt{2}, y/\sqrt{2} ),\] where \(C\) is a normalization constant that depends only on \(\boldsymbol{\textsl{k}}\).
In the following three propositions, we state several properties of the polynomials \(Q_n^{\boldsymbol{\textsl{k}}}(x,y)\).
Proposition 1. \({}\)
For \(n=0,1,\dots,k_l+1\) the polynomials \(Q_n^{\boldsymbol{\textsl{k}}}(x,y)\) satisfy \[Q^{\boldsymbol{\textsl{k}}}_n(x,y)=Q_n^{\boldsymbol{\textsl{k}}}(y,x)=Q^{\boldsymbol{\textsl{k}}}_n(-x,-y)=(-1)^{\kappa+n} Q^{\boldsymbol{\textsl{k}}}_n(-x,y),\] where \(\kappa:=k_1+k_2+\dots+k_l-l(l+1)/2\).
We have \[\label{eqn95Qk0} Q^{\boldsymbol{\textsl{k}}}_0(x,y)=(-2)^{l} \bigg[ \prod\limits_{j=1}^l k_j \bigg] \times H_{\boldsymbol{\textsl{k}}-1}(x) H_{\boldsymbol{\textsl{k}}-1}(y),\qquad{(2)}\] where \(\boldsymbol{\textsl{k}}-1:=(k_1-1,\dots,k_l-1)\).
For \(n=1,2,\dots,k_l+1\) the polynomials \(Q^{\boldsymbol{\textsl{k}}}_n(x,y)\) can be computed recursively via \[\label{Q95m95recursion} Q^{\boldsymbol{\textsl{k}}}_n(x,y)= {\mathbf{1}}_{\{n \notin \boldsymbol{\textsl{k}}\}} \frac{H_{\boldsymbol{\textsl{k}},n}(x) H_{\boldsymbol{\textsl{k}},n}(y)}{2^{n+l} n! \prod_{j=1}^l (n-k_j)}- \sum\limits_{m=1}^{n} \frac{H_{m}(x) H_{m}(y)}{2^{m} m!} Q^{\boldsymbol{\textsl{k}}}_{n-m}(x,y).\qquad{(3)}\]
The parity conditions and formula ?? were established in [1] for Krein-Adler sequences \(\boldsymbol{\textsl{k}}\).
Proposition 2. Let \(\boldsymbol{\textsl{k}}\in {\mathbb{N}}^l\) be a Krein-Adler sequence. Denote \[\begin{align} \Omega(x,y)&:=\frac{1}{y-x} \Big( \frac{H'_{\boldsymbol{\textsl{k}}}(y)}{H_{\boldsymbol{\textsl{k}}}(y)} - \frac{H'_{\boldsymbol{\textsl{k}}}(x)}{H_{\boldsymbol{\textsl{k}}}(x)} \Big), \\ \omega(x)&:=\lim_{y\to x} \Omega(x,y)= \frac{H_{\boldsymbol{\textsl{k}}}''(x)}{H_{\boldsymbol{\textsl{k}}}(x)}-\frac{H_{\boldsymbol{\textsl{k}}}'(x)^2}{H_{\boldsymbol{\textsl{k}}}(x)^2}. \end{align}\] For every \(x,y\in {\mathbb{C}}\), the following identities hold: \[\begin{align} \label{sum1} \sum\limits_{n=0}^{k_l+1} Q^{\boldsymbol{\textsl{k}}}_n(x,y)&= H_{\boldsymbol{\textsl{k}}}(x) H_{\boldsymbol{\textsl{k}}}(y),\\ \label{sum2} \sum\limits_{n=0}^{k_l+1} n Q_n^{\boldsymbol{\textsl{k}}}(x,y)&=- H_{\boldsymbol{\textsl{k}}}(x) H_{\boldsymbol{\textsl{k}}}(y)\big( \Omega(x,y)-l \big),\\ \label{sum3} \sum\limits_{n=0}^{k_l+1} n^2 Q_n^{\boldsymbol{\textsl{k}}}(x,y)&= H_{\boldsymbol{\textsl{k}}}(x) H_{\boldsymbol{\textsl{k}}}(y) \times \bigg((\Omega(x,y)-l)^2+\frac{\omega(x)+\omega(y)-2 \Omega(x,y)}{(y-x)^2} \bigg). \end{align}\] {#eq: sublabel=eq:sum1,eq:sum2,eq:sum3}
Formula ?? for general \(x\) and \(y\) and formula ?? for the special case \(x=y\) were first established in [1].
Proposition 3. Let \(k\in {\mathbb{N}}\) and \(\boldsymbol{\textsl{k}}=(k)\), so that \(l=1\). Then \[\label{Q95n95l611} Q_n^{\boldsymbol{\textsl{k}}}(x,y)={\mathbf{1}}_{\{n\geqslant 1\}} 2^{n-1} (k)_{n-1} H_{k+1-n}(x) H_{k+1-n}(y) - 2^{n+1} (k)_{n+1} H_{k-1-n}(x) H_{k-1-n}(y),\qquad{(4)}\] where \((r)_n:=r(r-1)\dots(r-n+1)\) denotes the falling factorial.
We also state several conjectures, which we have verified numerically for many cases of strictly increasing sequences \(\boldsymbol{\textsl{k}}\in {\mathbb{N}}^l\).
Conjecture 1: Identities ?? ?? hold for any strictly increasing sequence \(\boldsymbol{\textsl{k}}\in {\mathbb{N}}^l\), not only for Krein-Adler sequences.
Conjecture 2: For \(l\geqslant 2\) and any strictly increasing sequence \(\boldsymbol{\textsl{k}}\in {\mathbb{N}}^l\), the following identity holds: \[Q^{\boldsymbol{\textsl{k}}}_{k_l+1} (x,y)=2^{k_l+l-1} (k_l)! \bigg[ \prod\limits_{j=1}^{l-1} (k_l-k_j) \bigg] \times H_{\tilde{\boldsymbol{\textsl{k}}}}(x) H_{\tilde{\boldsymbol{\textsl{k}}}}(y),\] where \(\tilde{\boldsymbol{\textsl{k}}}=(k_1,k_2,\dots,k_{l-1})\).
Conjecture 3: When \(l=2\) and \(\boldsymbol{\textsl{k}}\) is a Krein-Adler sequence \(\boldsymbol{\textsl{k}}=[k,k+1]\) (with \(k \in {\mathbb{N}}\)), then \[\begin{align} Q^{\boldsymbol{\textsl{k}}}_{k+1}(x,y)&= 2^{k+1} k! k \Big( H_{k+1}(x) H_{k+1}(y) + 2 (k+1) \big[ H_{k+1}(x) H_{k-1}(y)\\ &\qquad \qquad +H_{k+1}(y) H_{k-1}(x) \big] + 4 k (k+1) H_{k-1}(x) H_{k-1}(y) \Big). \end{align}\]
Conjecture 4: The polynomials \(Q^{\boldsymbol{\textsl{k}}}_n(x,y)\) have integer coefficients. A stronger conjecture is that the coefficients of the polynomials \(Q^{\boldsymbol{\textsl{k}}}_n(x/2,y/2)\) are integers that are divisible by \[\prod\limits_{1\leqslant i < j \leqslant l} (k_i-k_j)^{2}.\]
Foata [13] gave a combinatorial proof of the classical Mehler formula. If true, Conjecture 4 may indicate that Theorem 1 could be proved by combinatorial methods and may suggest a potential combinatorial interpretation of the polynomials \(Q_n^{\boldsymbol{\textsl{k}}}(x,y)\).
We present the proofs of all our results in the next section. Our proof of Theorem 1 follows the same general strategy as the proof of Pupasov-Maksimov [1] in the case of Krein-Adler sequences. However, the Krein-Adler condition, which implies the orthogonality and completeness of \(f_{\boldsymbol{\textsl{k}},n}\) in \(L_2({\mathbb{R}}, {\textrm d}x)\), played an essential role in the proof in [1]. For the general case, we therefore had to modify the argument substantially and introduce several new ideas.
We denote by \(\psi_n(x):=e^{-x^2/2} H_n(x)\) the non-normalized Hermite functions. If \(n<0\) we set \(H_n(x)\equiv 0\) and \(\psi_n(x) \equiv 0\). For \(\boldsymbol{\textsl{j}}=(j_1,j_2,\dots,j_m) \in {\mathbb{Z}}_{\geqslant 0}^m\) we denote \(\psi_{\boldsymbol{\textsl{j}}}(x):=\textrm{Wr}[\psi_{j_1}(x), \dots,\psi_{j_m}(x)]\) and \[|\,\boldsymbol{\textsl{j}}\,|:=j_1+\dots+j_m.\] For a strictly increasing sequence \(\boldsymbol{\textsl{k}}\in {\mathbb{N}}^l\) and \(n \in {\mathbb{Z}}_{\geqslant 0}\) we denote \(\psi_{\boldsymbol{\textsl{k}},n}(x):=\psi_{\boldsymbol{\textsl{j}}}(x)\), where \(\boldsymbol{\textsl{j}}=(k_1,k_2,\dots,k_l,n)\). By the following well-known property of Wronskian determinants: \[\label{multiplication95identity} \textrm{Wr}[f(x) g_1(x),f(x) g_2(x) \dots,f(x) g_{m}(x)]=f(x)^m \textrm{Wr}[ g_1(x), g_2(x) \dots, g_{m}(x)],\tag{5}\] we have \[\label{psi95k95H95k} \psi_{\boldsymbol{\textsl{k}}}(x)=e^{-l x^2/2} H_{\boldsymbol{\textsl{k}}}(x), \;\;\; \psi_{\boldsymbol{\textsl{k}},n}(x)=e^{-(l+1) x^2/2} H_{\boldsymbol{\textsl{k}},n}(x).\tag{6}\] In what follows, we work with the bilinear generating functions of \(\psi_n\) and \(\psi_{\boldsymbol{\textsl{k}},n}\), defined by \[\label{def95mathcal95E} {\mathcal{E}}(z,x,y):=\sum\limits_{n\geqslant 0} z^n \frac{\psi_n(x) \psi_n(y)}{2^n n!}= \frac{1}{\sqrt{1-z^2}} \exp\bigg( \frac{4xyz-(1+z^2)(x^2+y^2)}{2(1-z^2)} \bigg),\tag{7}\] and \[\label{def95mathcal95Ek} {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y):=\sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} z^n \frac{\psi_{\boldsymbol{\textsl{k}},n}(x) \psi_{\boldsymbol{\textsl{k}},n}(y)}{2^{n+l} n! \prod\limits_{j=1}^l (n-k_j)}.\tag{8}\] From 6 , we have \[\label{mathcal95E95E} {\mathcal{E}}(z,x,y)=e^{-(x^2+y^2)/2} E(z,x,y), \;\;\; {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)=e^{-(l+1)(x^2+y^2)/2} E_{\boldsymbol{\textsl{k}}}(z,x,y).\tag{9}\]
Before proving Theorem 1, we need to establish several auxiliary results.
Lemma 1. \({}\)
There exists \(C=C(\boldsymbol{\textsl{k}})>0\) such that for all \(n\geqslant 1\) and \(x\in {\mathbb{C}}\) \[\label{psi95upper95bound} |\psi_{\boldsymbol{\textsl{k}},n}(x)|\leqslant C n^l \big | e^{-l x^2/2}\big | \times (1+|x|^{|\boldsymbol{\textsl{k}}|}) \times \sum\limits_{j=0}^l |\psi_{n-j}(x)|.\qquad{(5)}\]
The function \({\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)\) is an analytic function of three variables in the domain \(|z|<1\), \(x,y\in {\mathbb{C}}\).
The function \({\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)\) satisfies the following symmetry conditions: \[\label{E95k95symmetries} {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)={\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,y,x)={\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,-x,-y)=(-1)^{\kappa} {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(-z,-x,y), \;\;\; |z|<1, \; x,y\in {\mathbb{C}},\qquad{(6)}\] where \(\kappa=|\,\boldsymbol{\textsl{k}}\,|-l(l+1)/2\).
Proof. To prove (i), we start with the formula \[\psi_{\boldsymbol{\textsl{k}},n}(x)=e^{-(l+1)x^2/2} \textrm{Wr}[H_{k_1}(x),\dots,H_{k_l}(x),H_n(x)],\] and expand the Wronskian determinant along the last column, which contains the derivatives of \(H_n(x)\): \[\label{psi95k95n95sum95derivatives} \psi_{\boldsymbol{\textsl{k}},n}(x)=e^{-(l+1)x^2/2} \sum\limits_{j=0}^l p_{j}(x) \partial_x^j H_n(x),\tag{10}\] where \[p_{j}(x)=(-1)^{j+l} {\textrm{det}} \Big[ \partial_x^m H_{k_i}(x) \Big]_{\substack{1\leqslant i \leqslant l \\ 0\leqslant m \leqslant l \\ m\neq j}}.\] It is clear from the above determinant expression that \(p_{j}\) is a polynomial of degree not exceeding \(|\, \boldsymbol{\textsl{k}}\, |\). This, combined with the formula \[\partial_x^j H_m(x)=2^j m(m-1)\dots(m-j+1) H_{m-j}(x),\] implies the upper bound ?? .
The proof of item (ii) requires the following asymptotic result (see [14][Theorem 8.22.7]): as \(n\to +\infty\) \[\psi_n(x) = \frac{2^n}{\sqrt{\pi}} \Gamma((n+1)/2) \Big[ \cos(x \sqrt{2n+1} - n\pi/2)+O\big(n^{-1/2} \exp\big( |\textrm {Im}(x)|\sqrt{2n+1}\big)\big)\Big],\] uniformly in \(x\) on compact subsets of \({\mathbb{C}}\). We also need the fact that \[\frac{\Gamma((n+1)/2)^2}{n!}=2^{-n} \sqrt{\frac{2\pi}{n}} (1+o(1)), \;\;\; n\to +\infty,\] which follows from the Legendre duplication formula for the gamma function. Combining the above two facts with ?? , we conclude that for each compact subset \(A \subset {\mathbb{C}}\) there exists a constant \(C_1=C_1(A,\boldsymbol{\textsl{k}})>0\) such that the absolute value of each term in the infinite series 8 is bounded above by \[C_1 |z|^n n^{l-1/2} \exp\big(\sqrt{2n+1} (|\textrm {Im}(x)|+|\textrm {Im}(y)|)\big),\] for \(x,y \in A\). This implies that the series 8 converges uniformly on compact subsets of \(\{(z,x,y) \in {\mathbb{C}}^3 \; : \; |z|<1\}\) and defines an analytic function in this domain.
According to Lemma 2.1 and Lemma 3.6 in [15], the polynomial \(H_{\boldsymbol{\textsl{k}},n}(x)\) has degree \(\kappa+n\) and satisfies the parity condition \(H_{\boldsymbol{\textsl{k}},n}(-x)=(-1)^{\kappa+n} H_{\boldsymbol{\textsl{k}},n}(x)\). The symmetry conditions in item (iii) follow immediately from this result and formulas 6 , 7 and 8 . ◻
Lemma 2. Assume that \(U(x)\), \(g_1(x)\) and \(g_2(x)\) are smooth functions on \((-\infty, c)\) and that \[\label{f95i95equation} -g_i''(x)+U(x) g_i(x)=\lambda_i g_i(x), \;\;\; x \in (-\infty, c), \;\;\; i\in \{1,2\}.\qquad{(7)}\] Assume further that \(\lambda_1\neq \lambda_2\) and \(\textrm{Wr}[g_1(x),g_2(x)] \to 0\) as \(x\to -\infty\). Then, for any \(x\in (-\infty, c)\), we have \[\int_{-\infty}^x g_1(y)g_2(y) {\textrm d}y = \frac{1}{ \lambda_1-\lambda_2} \textrm{Wr}[g_1(x),g_2(x)].\]
Proof. The proof is based on the identity \[\frac{{\textrm d}}{{\textrm d}x}\textrm{Wr}[g_1(x),g_2(x)] = (\lambda_1-\lambda_2) g_1(x)g_2(x),\] which follows from ?? . ◻
Let \({\mathcal{L}}_{\boldsymbol{\textsl{k}},x}\) be the first-order differential operator acting on the \(x\)-variable, defined by \[{\mathcal{L}}_{\boldsymbol{\textsl{k}},x}h(x):=\textrm{Wr}[ \psi_{\boldsymbol{\textsl{k}}}(x), h(x)]=\psi_{\boldsymbol{\textsl{k}}}(x) h'(x) -\psi_{\boldsymbol{\textsl{k}}}'(x) h(x),\] where \(h\) is a smooth function on \({\mathbb{R}}\). For a strictly increasing sequence \(k=(k_1,\dots,k_l)\in {\mathbb{N}}^l\) we denote by \(\boldsymbol{\textsl{k}}\setminus k_j\) the sequence of length \(l-1\) obtained from \(\boldsymbol{\textsl{k}}\) by removing the term \(k_j\). We also denote \[\tilde{\boldsymbol{\textsl{k}}}:=\boldsymbol{\textsl{k}}\setminus k_l=(k_1,\dots,k_{l-1}).\] The Wronskian identity \[\label{Wronskian95identity} \textrm{Wr}[g_1,\dots,g_{j}]\times \textrm{Wr}[g_1,\dots,g_{j},h,f] =\textrm{Wr}\Big[\textrm{Wr}[g_1,\dots,g_{j},h],\textrm{Wr}[g_1,\dots,g_{j},f] \Big]\tag{11}\] implies the following key result \[\label{eqn95psi95kn95L95k} \psi_{\tilde{\boldsymbol{\textsl{k}}}}(x) \times \psi_{\boldsymbol{\textsl{k}},n}(x)= {\mathcal{L}}_{\boldsymbol{\textsl{k}},x} \psi_{\tilde{\boldsymbol{\textsl{k}}},n}(x).\tag{12}\]
Lemma 3. Let \(\boldsymbol{\textsl{k}}\in {\mathbb{N}}^l\) be a strictly increasing sequence. There exists \(c=c(\boldsymbol{\textsl{k}}) \in {\mathbb{R}}\) such that for \(|z|<1\), \(x\in {\mathbb{R}}\) and \(y \in (-\infty,c]\) we have \[\label{eqn95E95recursion} {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)=- \frac{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(y)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(x)} {\mathcal{L}}_{\boldsymbol{\textsl{k}},x} \int_{-\infty}^y {\mathcal{E}}_{\tilde{\boldsymbol{\textsl{k}}}}(z,x,u) \frac{\psi_{\boldsymbol{\textsl{k}}}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2} {\textrm d}u.\qquad{(8)}\]
Proof. We take \(c \in {\mathbb{R}}\) such that the polynomial \(H_{\tilde{\boldsymbol{\textsl{k}}}}(x)\) has no zeros on \((-\infty,c]\). This implies that \(\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)=\exp(-(l-1)u^2/2) H_{\tilde{\boldsymbol{\textsl{k}}}}(u)\) is also non-zero on \((-\infty,c]\). We recall that \(f_{\boldsymbol{\textsl{k}},n}(x)\) is defined in 2 and that these functions satisfy equation 3 .
To prove formula ?? , we first apply Lemma 2: \[\begin{align} \label{lemma395proof1} &\int_{-\infty}^y {\mathcal{E}}_{\tilde{\boldsymbol{\textsl{k}}}}(z,x,u) \frac{\psi_{\boldsymbol{\textsl{k}}}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2} {\textrm d}u= \sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \tilde{\boldsymbol{\textsl{k}}}} z^n \frac{\psi_{\tilde{\boldsymbol{\textsl{k}}},n}(x)}{2^{n+l-1} n! \prod\limits_{j=1}^{l-1} (n-k_j)} \int_{-\infty}^y f_{\tilde{\boldsymbol{\textsl{k}}},k_l}(u) f_{\tilde{\boldsymbol{\textsl{k}}},n}(u) {\textrm d}u\\ \nonumber &=\sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} z^n \frac{ \psi_{\tilde{\boldsymbol{\textsl{k}}},n}(x)}{2^{n+l-1} n! \prod\limits_{j=1}^{l-1} (n-k_j)} \frac{\textrm{Wr}[f_{\tilde{\boldsymbol{\textsl{k}}},k_l}(y),f_{\tilde{\boldsymbol{\textsl{k}}},n}(y)]}{2k_l-2n}+ \frac{z^{k_l} \psi_{\tilde{\boldsymbol{\textsl{k}}},k_l}(x)}{2^{k_l+l-1} k_l! \prod\limits_{j=1}^{l-1} (k_l-k_j)} \int_{-\infty}^y f_{\tilde{\boldsymbol{\textsl{k}}},k_l}(u)^2 {\textrm d}u. \end{align}\tag{13}\] From 12 we find \[\textrm{Wr}[f_{\tilde{\boldsymbol{\textsl{k}}},k_l}(y),f_{\tilde{\boldsymbol{\textsl{k}}},n}(y)]=\frac{1}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(y)^2} \textrm{Wr}[\psi_{\boldsymbol{\textsl{k}}}(y),\psi_{\tilde{\boldsymbol{\textsl{k}}},n}(y)]= \frac{\psi_{\boldsymbol{\textsl{k}},n}(y)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(y)},\] and therefore \[\int_{-\infty}^y {\mathcal{E}}_{\tilde{\boldsymbol{\textsl{k}}}}(z,x,u) \frac{\psi_{\boldsymbol{\textsl{k}}}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2} {\textrm d}u= -\frac{1}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(y)}\sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} z^n \frac{\psi_{\tilde{\boldsymbol{\textsl{k}}},n}(x) \psi_{\boldsymbol{\textsl{k}},n}(y)}{2^{n+l} n! \prod\limits_{j=1}^{l} (n-k_j)}+z^{k_l} \psi_{\boldsymbol{\textsl{k}}}(x) G(y),\] for some function \(G(y)\). Finally, applying identity 12 once again and noting that \({\mathcal{L}}_{\boldsymbol{\textsl{k}},x} \psi_{\boldsymbol{\textsl{k}}}(x)\equiv 0\), we obtain \[\begin{align} \label{lemma395proof2} - \frac{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(y)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(x)} {\mathcal{L}}_{\boldsymbol{\textsl{k}},x} \int_{-\infty}^y {\mathcal{E}}_{\tilde{\boldsymbol{\textsl{k}}}}(z,x,u) \frac{\psi_{\boldsymbol{\textsl{k}}}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2} {\textrm d}u&= \frac{1}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(x)}\sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} z^n \frac{{\mathcal{L}}_{\boldsymbol{\textsl{k}},x} \psi_{\tilde{\boldsymbol{\textsl{k}}},n}(x) \psi_{\boldsymbol{\textsl{k}},n}(y)}{2^{n+l} n! \prod\limits_{j=1}^{l} (n-k_j)} \\ \nonumber &= \sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} z^n \frac{\psi_{\boldsymbol{\textsl{k}},n}(x) \psi_{\boldsymbol{\textsl{k}},n}(y)}{2^{n+l} n! \prod\limits_{j=1}^{l} (n-k_j)}={\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y). \end{align}\tag{14}\]
In the above derivation of formula ?? , we interchanged summation and differentiation/integration. Let us first justify the interchange of summation and integration in the first step of 13 . Consider the function \(f_{\tilde{\boldsymbol{\textsl{k}}},k_l}(u) f_{\tilde{\boldsymbol{\textsl{k}}},n}(u)\), which we rewrite in the form \[f_{\tilde{\boldsymbol{\textsl{k}}},k_l}(u) f_{\tilde{\boldsymbol{\textsl{k}}},n}(u)=\Big[ e^{-u^2/2} \frac{H_{\boldsymbol{\textsl{k}}}(u)}{H_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2} \Big] \times \Big[ e^{(l-1)u^2/2} \psi_{\tilde{\boldsymbol{\textsl{k}}},n}(u)\Big].\] The above factorization follows from 2 and 6 . The function \(e^{-u^2/4} H_{\boldsymbol{\textsl{k}}}(u)/H_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2\) is bounded by some constant \(A=A(\boldsymbol{\textsl{k}})\) on \((-\infty, c]\) (since \(H_{\boldsymbol{\textsl{k}}}(u)\) and \(H_{\tilde{\boldsymbol{\textsl{k}}}}(u)\) are polynomials and \(H_{\tilde{\boldsymbol{\textsl{k}}}}(u)\) is non-zero on \((-\infty, c]\)). Thus, \(g(u):=e^{-u^2/2} H_{\boldsymbol{\textsl{k}}}(u)/H_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2\) satisfies \(|g(u)|<A e^{-u^2/4}\) for \(u \in (-\infty, c]\). From ?? we conclude that there exists a constant \(B=B(\tilde{\boldsymbol{\textsl{k}}})\) such that \[|e^{(l-1)u^2/2} \psi_{\tilde{\boldsymbol{\textsl{k}}},n}(u)|< B n^{l-1} (1+|u|^{M}) \sum\limits_{j=0}^{l-1} |\psi_{n-j}(u)|, \;\;\; u \in {\mathbb{R}},\] where \(M=|\tilde{\boldsymbol{\textsl{k}}}|\). Applying the Cauchy-Schwarz inequality and using the fact that \(\lVert \psi_n \rVert^2=\sqrt{\pi}2^n n!\), we obtain \[\begin{align} \int_{-\infty}^y | f_{\tilde{\boldsymbol{\textsl{k}}},k_l}(u) f_{\tilde{\boldsymbol{\textsl{k}}},n}(u) | {\textrm d}u &\leqslant A B n^{l-1} \sum\limits_{j=0}^{l-1} \int_{-\infty}^c (1+|u|^{M}) e^{-u^2/4} \; |\psi_{n-j}(u)|\, {\textrm d}u \\ &\leqslant A B n^{l-1} \Big[\int_{-\infty}^c (1+|u|^{M})^2 e^{-u^2/2}{\textrm d}u\Big]^{1/2} \times \sum\limits_{j=0}^{l-1} \lVert \psi_{n-j}\rVert \leqslant C l n^{l-1} \sqrt{2^n n!}, \end{align}\] for some \(C=C(\boldsymbol{\textsl{k}})>0\). Using the above estimate and the same argument as in the proof of Lemma 1, we conclude that the series \[\sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} r^n \frac{ |\psi_{\tilde{\boldsymbol{\textsl{k}}},n}(x)|}{2^{n+l-1} n! \prod\limits_{j=1}^{l-1} (n-k_j)} \int_{-\infty}^y | f_{\tilde{\boldsymbol{\textsl{k}}},k_l}(u) f_{\tilde{\boldsymbol{\textsl{k}}},n}(u) | {\textrm d}u\] converges for all \(r\in (0,1)\), \(x\in {\mathbb{C}}\) and \(y \in (-\infty, c]\). Thus, we can apply Fubini’s Theorem and interchange summation and integration in 13 .
To justify the interchange of summation and differentiation in the first step of 14 , we use the same argument as in the proof of Lemma 1 and conclude that the series \[\sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} z^n \frac{\psi_{\tilde{\boldsymbol{\textsl{k}}},n}(x) \psi_{\boldsymbol{\textsl{k}},n}(y)}{2^{n+l} n! \prod\limits_{j=1}^{l} (n-k_j)}\] converges uniformly on compact subsets of \(\{(z,x,y) \in {\mathbb{C}}^3 \; : \; |z|<1\}\). Hence, it defines an analytic function on this set, and the infinite series can be differentiated term-by-term. ◻
For \(n\in {\mathbb{N}}\), \(|z|<1\), \(x\in {\mathbb{C}}\) and \(y\in {\mathbb{R}}\) we denote \[\label{def95In} I_n(z,x,y):=\int_{-\infty}^y {\mathcal{E}}(z,x,u) \psi_n(u) {\textrm d}u.\tag{15}\]
Lemma 4. The function \(I_n(z,x,y)\) extends to an analytic function in the domain \(\{(z,x,y) \in {\mathbb{C}}^3 \; : \; |z|<1\}\) and it satisfies \[\label{eqn95I95n95main} I_n(z,x,y)= z^n \psi_n(x) \frac{\sqrt{\pi}}{2} \big( 1+ {\textrm{erf}}(w) \big)-e^{-x^2/2-w^2} \sum\limits_{i=1}^n \binom{n}{i} \zeta^{i} z^{n-i} H_{i-1}(w) H_{n-i}(x),\qquad{(9)}\] where \(\zeta:=\sqrt{1-z^2}\) and \(w:=(y-xz)/\zeta\).
Proof. We first observe that \({\mathcal{E}}(z,x,y)\) defined in 7 can be written in the form \[{\mathcal{E}}(z,x,y)=\zeta^{-1} e^{y^2/2-x^2/2-(y-xz)^2/\zeta^2}.\] Assume that \(z\in (-1,1)\) and \(x,y\in {\mathbb{R}}\). Using the above expression for \({\mathcal{E}}(z,x,y)\), we rewrite 15 in the equivalent form \[I_n(z,x,y)=\zeta^{-1} e^{-x^2/2} \int_{-\infty}^y e^{-(u-xz)^2/\zeta^2} H_n(u) {\textrm d}u.\] We change the variable of integration to \(v=(u-xz)/\zeta\) and obtain \[I_n(z,x,y)=e^{-x^2/2} \int_{-\infty}^w e^{-v^2} H_n(\zeta v + x z) {\textrm d}v.\] Noting that \(\zeta^2+z^2=1\), we apply the addition formula for Hermite polynomials (see [2]) \[H_n(\zeta v+xz)=\sum\limits_{i=0}^n \binom{n}{i} \zeta^i z^{n-i} H_i(v) H_{n-i}(x).\] The desired result ?? follows from the above two equations and the following formulas: \[\label{two95integrals} \int_{-\infty}^w e^{-v^2} {\textrm d}v= \frac{\sqrt{\pi}}{2} \big( 1+ {\textrm{erf}}(w) \big), \;\;\; \int_{-\infty}^w e^{-v^2} H_i(v) {\textrm d}v=-e^{-w^2} H_{i-1}(w),\tag{16}\] which can be found in [2] and [16].
Now that we have established ?? for \(z \in (-1,1)\), \(x,y\in {\mathbb{R}}\), we note that the right-hand side in ?? is an analytic function in \(\{(z,x,y) \in {\mathbb{C}}^3 \; : \; |z|<1\}\). Thus, \(I_n(z,x,y)\) can be extended in this domain by analytic continuation. ◻
We introduce the \(l\)-th order linear differential operator \({\mathcal{M}}_{\boldsymbol{\textsl{k}},x}\), acting on the \(x\)-variable, by \[\label{def95M95kx} {\mathcal{M}}_{\boldsymbol{\textsl{k}},x} h(x):=\textrm{Wr}[ \psi_{k_1},\dots,\psi_{k_l},h](x),\tag{17}\] where \(h\) is an entire function. The following result was the starting point for the derivation of the Mehler formula for exceptional Hermite polynomials in [1]; see also [17], where this result was first obtained for Krein-Adler sequences \(\boldsymbol{\textsl{k}}\).
Lemma 5. For \(|z|<1\), \(x,y\in {\mathbb{C}}\) we have \[\label{eqn95PM} {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)= {\mathcal{M}}_{\boldsymbol{\textsl{k}},x} \sum\limits_{j=1}^l (-1)^{j+l-1} \psi_{\boldsymbol{\textsl{k}}\setminus k_j}(y) I_{k_j}(z,x,y).\qquad{(10)}\]
Proof. The proof proceeds by induction on \(l\). The case \(l=1\) and \(\boldsymbol{\textsl{k}}=(k)\) is covered by Lemma 3, followed by analytic continuation in \(x\) and \(y\). Assume now that \(l\geqslant 2\) and that for \(\tilde{\boldsymbol{\textsl{k}}}=(k_1,k_2,\dots,k_{l-1})\) we have \[{\mathcal{E}}_{\tilde{\boldsymbol{\textsl{k}}}}(z,x,y)= {\mathcal{M}}_{\tilde{\boldsymbol{\textsl{k}}},x} \sum\limits_{j=1}^{l-1} (-1)^{j+l-2} \psi_{\tilde{\boldsymbol{\textsl{k}}}\setminus k_j}(y) I_{k_j}(z,x,y),\] with the above identity valid for all \(z \in (-1,1)\) and \(x,y\in {\mathbb{R}}\). We substitute the above expression for \({\mathcal{E}}_{\tilde{\boldsymbol{\textsl{k}}}}(z,x,y)\) into the right-hand side of ?? and conclude that there exists \(c \in {\mathbb{R}}\) such that for \(z\in (-1,1)\), \(x\in {\mathbb{R}}\) and \(y \in (-\infty,c]\) we have \[\label{lemma595proof1} {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)= \frac{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(y)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(x)} {\mathcal{L}}_{\boldsymbol{\textsl{k}},x}{\mathcal{M}}_{\tilde{\boldsymbol{\textsl{k}}},x} \sum\limits_{j=1}^{l-1} (-1)^{j+l-1} \int_{-\infty}^y \psi_{\tilde{\boldsymbol{\textsl{k}}}\setminus k_j}(u) I_{k_j}(z,x,u) \frac{\psi_{\boldsymbol{\textsl{k}}}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2} {\textrm d}u.\tag{18}\] Using the Wronskian identity 11 , we check that \[\frac{{\textrm d}}{{\textrm d}u} \frac{\psi_{\boldsymbol{\textsl{k}}\setminus k_j}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2}=\psi_{\tilde{\boldsymbol{\textsl{k}}}\setminus k_j}(u) \frac{\psi_{\boldsymbol{\textsl{k}}}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2}.\] Note that \(\psi_{\boldsymbol{\textsl{k}}\setminus k_j}(u)/\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)\) is a rational function of \(u\) (this follows from 5 ). Due to our assumption \(z\in (-1,1)\), we have \(\zeta>0\), thus from ?? we see that the function \(u\mapsto I_{k_j}(z,x,u)\) decays exponentially fast as \(u\to -\infty\) (this fact also follows from the integral definition 15 ). Thus, \[\lim\limits_{u\to -\infty} I_{k_j}(z,x,u)\frac{\psi_{\boldsymbol{\textsl{k}}\setminus k_j}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)}=0.\] Using the above two results and integrating by parts, we obtain \[\label{lemma595proof2} \int_{-\infty}^y \psi_{\tilde{\boldsymbol{\textsl{k}}}\setminus k_j}(u) I_{k_j}(z,x,u) \frac{\psi_{\boldsymbol{\textsl{k}}}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)^2} {\textrm d}u= I_{k_j}(z,x,y) \frac{\psi_{\boldsymbol{\textsl{k}}\setminus k_j}(y)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(y)}-\int_{-\infty}^y {\mathcal{E}}(z,x,u) \psi_{k_j}(u) \frac{\psi_{\boldsymbol{\textsl{k}}\setminus k_j}(u)}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(u)} {\textrm d}u.\tag{19}\]
Next, we use the identity \[0=\det \begin{bmatrix} \psi_{k_1}(u) & \psi_{k_2}(u) & \cdots & \psi_{k_{l}}(u)\\ \psi_{k_1}(u) & \psi_{k_2}(u) & \cdots & \psi_{k_{l}}(u)\\ \psi_{k_1}'(u) & \psi_{k_2}'(u) & \cdots & \psi_{k_{l}}'(u)\\ \vdots & \vdots & \ddots & \vdots\\ \psi_{k_1}^{(l-2)}(u) & \psi_{k_2}^{(l-2)}(u) & \cdots & \psi_{k_{l}}^{(l-2)}(u) \end{bmatrix}=\sum\limits_{j=1}^l (-1)^{j-1} \psi_{k_j}(u)\times \psi_{\boldsymbol{\textsl{k}}\setminus k_j}(u),\] where we have expanded the above determinant along the first row. Therefore, \[\sum\limits_{j=1}^{l-1} (-1)^{j-1} \psi_{k_j}(u) \psi_{\boldsymbol{\textsl{k}}\setminus k_j}(u) = (-1)^{l} \psi_{k_l}(u) \psi_{\tilde{\boldsymbol{\textsl{k}}}}(u).\] Combining the above result with 18 and 19 we arrive at \[{\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)= \frac{1}{\psi_{\tilde{\boldsymbol{\textsl{k}}}}(x)} {\mathcal{L}}_{\boldsymbol{\textsl{k}},x}{\mathcal{M}}_{\tilde{\boldsymbol{\textsl{k}}},x} \sum\limits_{j=1}^{l} (-1)^{j+l-1} \psi_{\boldsymbol{\textsl{k}}\setminus k_j}(y) I_{k_j}(z,x,y).\] From the Wronskian identity 11 it follows that \[{\mathcal{L}}_{\boldsymbol{\textsl{k}},x}{\mathcal{M}}_{\tilde{\boldsymbol{\textsl{k}}},x} h(x) = \psi_{\tilde{\boldsymbol{\textsl{k}}}}(x) {\mathcal{M}}_{\boldsymbol{\textsl{k}},x}h(x),\] for any entire function \(h\). Thus, we have proved that for \(z\in (-1,1)\), \(x\in {\mathbb{R}}\) and \(y\in (-\infty, c]\) we have \[{\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)= {\mathcal{M}}_{\boldsymbol{\textsl{k}},x} \sum\limits_{j=1}^{l} (-1)^{j+l-1} \psi_{\boldsymbol{\textsl{k}}\setminus k_j}(y) I_{k_j}(z,x,y).\] By analytic continuation, the above identity holds everywhere in the domain \(\{(z,x,y)\in {\mathbb{C}}^3 \; : \; |z|<1\}\). This completes the induction step. ◻
Corollary 1. For \(l\geqslant 2\), \(|z|<1\), \(x\in {\mathbb{C}}\) and \(y\in {\mathbb{R}}\) we have \[\label{eqn95PM2} {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)= {\mathcal{M}}_{\boldsymbol{\textsl{k}},x} \int_{-\infty}^y {\mathcal{E}}(z,x,u) e^{-u^2/2} P_{\boldsymbol{\textsl{k}}}(u,y) {\textrm d}u,\qquad{(11)}\] where \(P_{\boldsymbol{\textsl{k}}}(u,y)\) is a polynomial in \(u\) that satisfies \(\partial_u^j P_{\boldsymbol{\textsl{k}}}(u,y) |_{u=y}=0\) for all \(j=0,1,\dots,l-2\).
Proof. Formula ?? is equivalent to ?? with \(P_{\boldsymbol{\textsl{k}}}(u,y)\) defined by \[e^{-u^2/2} P_{\boldsymbol{\textsl{k}}}(u,y)=\sum\limits_{j=1}^l (-1)^{j+l-1} \psi_{k_j}(u) \psi_{\boldsymbol{\textsl{k}}\setminus k_j}(y)=(-1)^l \det \begin{bmatrix} \psi_{k_1}(u) & \psi_{k_2}(u) & \cdots & \psi_{k_{l}}(u)\\ \psi_{k_1}(y) & \psi_{k_2}(y) & \cdots & \psi_{k_{l}}(y)\\ \psi_{k_1}'(y) & \psi_{k_2}'(y) & \cdots & \psi_{k_{l}}'(y)\\ \vdots & \vdots & \ddots & \vdots\\ \psi_{k_1}^{(l-2)}(y) & \psi_{k_2}^{(l-2)}(y) & \cdots & \psi_{k_{l}}^{(l-2)}(y) \end{bmatrix}.\] The above equation can be written in the form \[P_{\boldsymbol{\textsl{k}}}(u,y)=e^{-(l-1)y^2/2} (-1)^l \det\begin{bmatrix} H_{k_1}(u) & H_{k_2}(u) & \cdots & H_{k_{l}}(u)\\ H_{k_1}(y) & H_{k_2}(y) & \cdots & H_{k_{l}}(y)\\ H_{k_1}'(y) & H_{k_2}'(y) & \cdots & H_{k_{l}}'(y)\\ \vdots & \vdots & \ddots & \vdots\\ H_{k_1}^{(l-2)}(y) & H_{k_2}^{(l-2)}(y) & \cdots & H_{k_{l}}^{(l-2)}(y) \end{bmatrix},\] which immediately implies \(\partial_u^j P_{\boldsymbol{\textsl{k}}}(u,y) |_{u=y}=0\) for all \(j=0,1,\dots,l-2\). ◻
Lemma 6. Let \(m,n \in {\mathbb{Z}}_{\geqslant 0}\) and let \(x, y\) be fixed real numbers satisfying \(x>y\). Let \(P(u)\) be a polynomial such that \(P^{(i)}(y)=0\) for \(i=0,1,\dots,n\). If \(m\leqslant n+1\) then \[\label{eqn95integral95limit} \lim\limits_{z\to 1-} \frac{\partial_x^{m} \int_{-\infty}^y {\mathcal{E}}(z,x,u) e^{-u^2/2} P(u) {\textrm d}u}{{\mathcal{E}}(z,x,y)} =0.\qquad{(12)}\] If \(m=n+2\), the above limit exists and is finite.
Proof. First, we check that formula 7 can be rewritten in the equivalent form \[\label{E95factorization} {\mathcal{E}}(z,x,y) =e^{(x^2+y^2)/2 - 2xy/(1+z)} \times \zeta^{-1} e^{-(y-x)^2/\zeta^2 },\tag{20}\] where \(\zeta=\zeta(z)=\sqrt{1-z^2}\). For the rest of this proof, we assume that \(1/\sqrt{2}<z<1\), so that \(0<\zeta<1/\sqrt{2}\).
We claim that, for any \(m\in {\mathbb{Z}}_{\geqslant 0}\) and any fixed \(x,y \in {\mathbb{R}}\) such that \(x>y\), we have \[\label{eqn95q95integral} \zeta^{-2} e^{(y-x)^2/\zeta^2} \int_{-\infty}^y e^{-2xu/(1+z)-(u-x)^2/\zeta^2} u^m {\textrm d}u \to \frac{e^{-x y}y^m}{2(x-y)},\tag{21}\] as \(z\to 1-\) (note that \(\zeta\to 0+\) as \(z\to 1-\)). Indeed, changing the variable of integration \(u=y-\zeta^2 v\), we obtain \[\zeta^{-2} e^{(y-x)^2/\zeta^2} \int_{-\infty}^y e^{-2xu/(1+z)-(u-x)^2/\zeta^2} u^m {\textrm d}u = e^{-2x y/(1+z)} \int_0^{\infty} e^{-2(x-y)v -\zeta^2 (v^2-2xv/(1+z))} (y-\zeta^2 v)^m {\textrm d}v,\] and the limit in 21 follows by the dominated convergence theorem.
Next, let \(q(z,u)\) be a polynomial in \(u\) of the form \[\label{form95of95q} q(z,u)=\sum\limits_{i=0}^M a_i(z) u^i,\tag{22}\] where \(M\) does not depend on \(z\) and the coefficients \(a_i(z)\) are continuous in some neighborhood of \(z=1\). For \(z\in (0,1)\) and \(j \in {\mathbb{Z}}_{\geqslant 0}\) we consider \[\label{def95Q95j} {\mathcal{Q}}_{j}(z,x,y):=\zeta^{-1} \int_{-\infty}^y \Big[\partial_u^j e^{-(u-x)^2/\zeta^2} \Big] \times e^{-2xu/(1+z)} q(z,u) (u-y)^{n+1} {\textrm d}u.\tag{23}\] From 20 and 21 we conclude that \[\frac{{\mathcal{Q}}_{0}(z,x,y)}{{\mathcal{E}}(z,x,y)} \to 0\] as \(z \to 1-\). When \(0< j \leqslant n+1\), we integrate by parts \(j\) times and obtain \[\begin{align} \label{Q95j95integration95by95parts} {\mathcal{Q}}_j(z,x,y)&= \zeta^{-1} \sum\limits_{i=0}^{j-1} (-1)^i \Big[\partial_u^{j-1-i} e^{-(u-x)^2/\zeta^2} \Big]_{u=y} \times \Big[\partial_u^i \Big( e^{-2xu/(1+z)} q(z,u) (u-y)^{n+1}\Big)\Big]_{u=y} \\\nonumber &+(-1)^j \zeta^{-1} \int_{-\infty}^y e^{-(u-x)^2/\zeta^2} \partial_u^j \Big( e^{-2xu/(1+z)} q(z,u) (u-y)^{n+1}\Big) {\textrm d}u. \end{align}\tag{24}\] In the finite sum in the above formula, we have \(i\leqslant j-1 \leqslant n\), and therefore \[\Big[\partial_u^i e^{-2xu/(1+z)} q(z,u) (u-y)^{n+1}\Big]_{u=y}=0.\] Thus, the finite sum in 24 is equal to zero, and it follows from 20 and 21 that, for \(0<j\leqslant n+1\), \[\frac{{\mathcal{Q}}_{j}(z,x,y)}{{\mathcal{E}}(z,x,y)} \to 0,\] as \(z \to 1-\).
When \(j=n+2\), the preceding argument shows that the terms with \(0\leqslant i \leqslant j-2\) in the finite sum in 24 are all equal to zero, and when \(i=j-1=n+1\) we have \[\Big[\partial_u^{n+1} \Big( e^{-2xu/(1+z)} q(z,u) (u-y)^{n+1}\Big)\Big]_{u=y}= e^{-2xy/(1+z)} q(z,y) (n+1)!.\] Thus, it follows from 20 , 21 and 24 that, as \(z\to 1-\), \[\frac{{\mathcal{Q}}_{n+2}(z,x,y)}{{\mathcal{E}}(z,x,y)} \to (-1)^{n+1} e^{-(x^2+y^2)/2} q(1,y) (n+1)!.\]
Now we are ready to complete the proof of Lemma 6. For any fixed \(z\in (1/\sqrt{2},1)\), we interchange the integration and differentiation in the numerator in ?? , use the factorization 20 , and obtain \[\label{numerator95as95sum95of95Q} \partial_x^{m} \int_{-\infty}^y {\mathcal{E}}(z,x,u) e^{-u^2/2} P(u) {\textrm d}u = \sum\limits_{j=0}^m (-1)^j \binom{m}{j} \widetilde{\mathcal{Q}}_{j}(z,x,y),\tag{25}\] where we have defined \[\label{def95tilde95Q} \widetilde{\mathcal{Q}}_{j}(z,x,y):=\zeta^{-1} \int_{-\infty}^y \Big[ \partial_u^{j} e^{-(u-x)^2/\zeta^2} \Big] \times \Big[ \partial_x^{m-j} e^{x^2/2-2xu/(1+z)} \Big] \times P(u) {\textrm d}u.\tag{26}\] When deriving the above formula, we used the fact that \[\partial_x^{j} e^{-(u-x)^2/\zeta^2}=(-1)^j \partial_u^{j} e^{-(u-x)^2/\zeta^2}.\] The interchange of integration and differentiation is justified since, for \(z\in (1/\sqrt{2},1)\), we have \(1/\zeta^2>2\); thus, for all \(u\in {\mathbb{R}}\) the integrand in 26 is dominated by \(C\exp(- u^2)\) for some \(C>0\), uniformly in \(x\) on compact subsets of \({\mathbb{R}}\). We note that the conditions \(P^{(i)}(y)=0\) for \(i=0,1,\dots,n\) imply that \(P(u)=q_1(u) (u-y)^{n+1}\) for some polynomial \(q_1\), which in turn implies \[\Big[ \partial_x^{m-j} e^{x^2/2-2xu/(1+z)} \Big] \times P(u)= e^{x^2/2-2xu/(1+z)} q_2(x,u/(1+z),u) (u-y)^{n+1},\] where \(q_2(x_1,x_2,x_3)\) is some polynomial in three variables. The important conclusion is that \(q_3(z,u)=q_2(x,u/(1+z),u)\) is of the form 22 (recall that everywhere in this proof \(x\) is fixed, so dependence of the coefficients of the polynomial \(q_3(z,u)\) on \(x\) is not a concern). Thus, \(\widetilde{\mathcal{Q}}_{j}(z,x,y)\) is an integral of the form 23 , and using the results about \({\mathcal{Q}}_j(z,x,y)\) that we established above, we conclude that the limit \[\lim\limits_{z\to 1-} \frac{\widetilde{\mathcal{Q}}_{j}(z,x,y)}{{\mathcal{E}}(z,x,y)}\] exists and is equal to zero for all \(0\leqslant j \leqslant n+1\), while for \(j=n+2\) the limit exists and is finite. The desired result now follows from 25 . ◻
Corollary 2. Let \(x\) and \(y\) be fixed real numbers such that \(x>y\). Then \[\lim_{z\to 1-} \frac{{\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)}{{\mathcal{E}}(z,x,y)}\] exists and is finite.
Proof. Using the same reasoning as in the derivation of 10 , we conclude that \[\label{M95kx95h} {\mathcal{M}}_{\boldsymbol{\textsl{k}},x} h(x)=e^{-l x^2/2} \sum\limits_{j=0}^l q_{j}(x) \partial_x^j h(x),\tag{27}\] for some real polynomials \(q_j(x)\), which may also depend on \(\boldsymbol{\textsl{k}}\). The desired result now follows from Corollary 1 and Lemma 6. ◻
Proof of Theorem 1: We denote by \({\mathbb{R}}[x]\) the set of polynomials in the \(x\)-variable with real coefficients, and by \({\mathbb{R}}_{n}[z,x,y]\) the set of real polynomials in the three variables \(z,x,y\) whose degree in \(z\) is not greater than \(n\). An expression such as \[g(z,x,y) \in H(z,x,y) \sum\limits_{i=1}^n {\mathbb{R}}[x] h_{i}(z,x,y)\] should be interpreted as stating that there exist polynomials \(P_i \in {\mathbb{R}}[x]\) such that \[g(z,x,y) = H(z,x,y) \sum\limits_{i=1}^n P_i(x) h_{i}(z,x,y)\] for all \(z,x,y\) in the domain of \(g\), \(H\) and \(h_i\). For example, with this notation we can express 27 in the form \[\label{item95i} {\mathcal{M}}_{\boldsymbol{\textsl{k}},x} h(x) \in e^{-l x^2/2}\sum\limits_{j=0}^l {\mathbb{R}}[x] \, \partial_x^j h(x).\tag{28}\]
As before, we denote \(\zeta:=\sqrt{1-z^2}\) and \(w:=(y-xz)/\zeta\). We collect some preliminary facts.
For all \(n\in {\mathbb{Z}}_{\geqslant 0}\) \[\label{item95ii} H_n(w)=n! \sum\limits_{j=0}^{\lfloor n/2 \rfloor} \frac{(-1)^j (2w)^{n-2j}}{j! (n-2j)!} \in \zeta^{-n} {\mathbb{R}}_{n}[z,x,y].\tag{29}\]
For all \(n\in {\mathbb{N}}\) \[\label{item95iii} \partial_x^n \, {\textrm{erf}}(w)=-(z/\zeta)^n \frac{2}{\sqrt{\pi}} e^{-w^2} H_{n-1}(w) \in e^{-w^2} \zeta^{1-2n} {\mathbb{R}}_{2n-1}[z,x,y].\tag{30}\]
For all \(n,m \in {\mathbb{Z}}_{\geqslant 0}\) \[\label{item95iv} \partial_x^m \Big[ e^{-w^2} H_n(w)\Big]=(z/\zeta)^m e^{-w^2} H_{n+m}(w).\tag{31}\]
The expression for \(H_n(w)\) in item (i) can be found in [2]. The results in items (ii) and (iii) follow from 16 and the fact that \(\partial_x = (-z/\zeta) \partial_w\).
We denote \[\begin{align} I^{(1)}_n(z,x,y)&:= z^n \psi_n(x) \frac{\sqrt{\pi}}{2} \big( 1+ {\textrm{erf}}(w) \big), \\ I_n^{(2)}(z,x,y)&:=-z^n e^{-w^2} \sum\limits_{i=1}^n \binom{n}{i} (z/\zeta)^{-i} H_{i-1}(w) \psi_{n-i}(x). \end{align}\] By ?? , we have \[{\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y) \in e^{-(l-1) y^2/2} \sum\limits_{j=1}^l {\mathbb{R}}[y] \Big[ {\mathcal{M}}_{\boldsymbol{\textsl{k}},x} I^{(1)}_{k_j}(z,x,y)+{\mathcal{M}}_{\boldsymbol{\textsl{k}},x} I^{(2)}_{k_j}(z,x,y)\Big].\] Using 28 we obtain \[{\mathcal{M}}_{\boldsymbol{\textsl{k}},x} I^{(1)}_{k_j}(z,x,y) \in z^{k_j} \frac{\sqrt{\pi}}{2} \big( 1+ {\textrm{erf}}(w) \big) {\mathcal{M}}_{\boldsymbol{\textsl{k}},x} \psi_{k_j}(x)+ z^{k_j} e^{-lx^2/2} \sum\limits_{\substack{n\geqslant 0, m\geqslant 1 \\ m+n\leqslant l}} {\mathbb{R}}[x] \partial_x^n \psi_{k_j}(x)\partial_x^m {\textrm{erf}}(w).\] By definition of \({\mathcal{M}}_{\boldsymbol{\textsl{k}},x}\) (see 17 ), we have \({\mathcal{M}}_{\boldsymbol{\textsl{k}},x} \psi_{k_j}(x)\equiv 0\). We use 29 , 30 , and the fact that \(\partial_x^n \psi_{k_j}(x)\in \exp(-x^2/2) {\mathbb{R}}[x]\) to conclude that \[{\mathcal{M}}_{\boldsymbol{\textsl{k}},x} I^{(1)}_{k_j}(z,x,y) \in z^{k_j} e^{-(l+1)x^2/2-w^2} \sum\limits_{m=1}^l {\mathbb{R}}[x] (z/\zeta)^m H_{m-1}(w)\subset e^{-(l+1)x^2/2-w^2} \zeta^{1-2l} {\mathbb{R}}_{2l-1+k_j}[z,x,y].\] In the last step, we used the fact that for all \(m=1,2,\dots,l\) \[z^{k_j+m} \zeta^{-m} H_{m-1}(w) \in z^{k_j+m} \zeta^{1-2m} {\mathbb{R}}_{m-1}[z,x,y] \subset \zeta^{1-2l} {\mathbb{R}}_{2l-1+k_j}[z,x,y],\] which is true since \(\zeta^{2l-2m}=(1-z^2)^{l-m}\). Similarly, using 29 and 31 , we conclude that \[{\mathcal{M}}_{\boldsymbol{\textsl{k}},x} I^{(2)}_{k_j}(z,x,y) \in z^{k_j} e^{-(l+1)x^2/2-w^2} \sum\limits_{i=1}^{k_j} \sum\limits_{m=0}^l {\mathbb{R}}[x] (z/\zeta)^{m-i} H_{m+i-1}(w) \subset e^{-(l+1)x^2/2-w^2} \zeta^{1-2l} {\mathbb{R}}_{2l-1+k_j}[z,x,y].\] Thus, \[{\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y) \in e^{-(l+1)x^2/2-(l-1)y^2/2-w^2} \zeta^{1-2l} {\mathbb{R}}_{2l-1+k_l}[z,x,y].\] Using the expression 20 and the above two equations we arrive at \[\frac{{\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)}{{\mathcal{E}}(z,x,y)} \in e^{-l(x^2+y^2)/2} (1-z^2)^{1-l} \, {\mathbb{R}}_{2l-1+k_l}[z,x,y].\] Thus, we have established that \[\label{R95and95P} R(z,x,y):=e^{l(x^2+y^2)/2} \frac{{\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)}{{\mathcal{E}}(z,x,y)}=\frac{P(z,x,y)}{(1-z^2)^{l-1}},\tag{32}\] where \(P(z,x,y)\) is a real polynomial in three variables, whose degree in the \(z\)-variable is at most \(2l-1+k_l\). We see that \(R(z,x,y)\) is a rational function of \(z\) with possible poles only at \(z=\pm 1\). Now assume that \(x, y\) are real and \(x>y\). Letting \(z\to 1-\) and applying Corollary 2, we conclude that \(R(z,x,y)\) has a finite limit at \(z=1\). This implies that, when \(x, y\) are real and \(x>y\), the polynomial \(P(z,x,y)\) is divisible by \((1-z)^{l-1}\). This divisibility property holds for all \(x,y\in {\mathbb{C}}\). Indeed, consider the Taylor coefficients of the numerator \(P(z,x,y)\) at \(z=1\). These Taylor coefficients are polynomials in \(x,y\) and the fact that \(P(z,x,y)\) is divisible by \((1-z)^{l-1}\) implies that the first \(l-1\) coefficients are equal to zero. It is clear that if a polynomial in two variables \(x,y\) is zero for all real values \(x>y\), then it must be equal to zero for all \(x,y\in {\mathbb{C}}\).
By the symmetry conditions ?? , the rational function \(R\) satisfies \(R(z,x,y)=(-1)^{\kappa} R(-z,-x,y)\), which implies that it also has a finite limit as \(z\to -1\). Thus, the polynomial \(P(z,x,y)\) is divisible by \((1-z^2)^{l-1}\). Since \((1-z^2)^{l-1}\) has degree \(2l-2\), and the numerator \(P(z,x,y)\) has degree in \(z\) not greater than \(2l-1+k_l\), we conclude that the function \(R(z,x,y)\) must be a polynomial in the three variables \(z,x,y\), with degree in \(z\) not greater than \(k_l+1\).
Thus, we have established that there exists a polynomial \(R(z,x,y)\), with degree in \(z\) not greater than \(k_l+1\), such that \[{\mathcal{E}}_{\boldsymbol{\textsl{k}}}(z,x,y)=e^{-l(x^2+y^2)/2}\, {\mathcal{E}}(z,x,y) R(z,x,y),\] which is equivalent to the desired result ?? by 9 . \(\sqcap\kern-8.0pt\sqcup\)
Remark 2: One possible way to prove Conjecture 2 is to identify the coefficient of the highest power of \(z\) in the polynomial \(P(z,x,y)\) in 32 (namely, the coefficient of \(z^{2l-1+k_l}\)). Indeed, after dividing \(P(z,x,y)\) by \((1-z^2)^{l-1}\), this would give \(Q^{\boldsymbol{\textsl{k}}}_{k_l+1}(x,y)\) – the coefficient of \(z^{k_l+1}\) in the polynomial function \(R(z,x,y)\). Thus, in order to prove Conjecture 2, one
would need to follow the steps of the proof of Theorem 1 and keep track of the coefficient of the highest power of \(z\) in every term. However,
the resulting expressions become complicated very quickly, and we were unable to complete this computation.
Proof of Proposition 1: We rewrite ?? in the form \[\label{generating95functions} \sum\limits_{n \in {\mathbb{Z}}_{\geqslant 0} \setminus \boldsymbol{\textsl{k}}} z^n \frac{H_{\boldsymbol{\textsl{k}},n}(x) H_{\boldsymbol{\textsl{k}},n}(y)}{2^{n+l} n! \prod\limits_{j=1}^l (n-k_j)}= \sum\limits_{n\geqslant 0} z^n \frac{H_n(x) H_n(y)}{2^n n!} \sum_{m=0}^{k_l+1} z^m Q_m^{\boldsymbol{\textsl{k}}}(x,y).\tag{33}\] Comparing the coefficient of \(z^0\), we obtain \[Q_0^{\boldsymbol{\textsl{k}}}(x,y)= \frac{H_{\boldsymbol{\textsl{k}},0}(x) H_{\boldsymbol{\textsl{k}},0}(y)}{(-2)^{l} \prod\limits_{j=1}^l k_j}.\] Since \(H_{\boldsymbol{\textsl{k}},0}(x)=(-2)^l \Big[\prod_{j=1}^l k_j\Big] H_{\boldsymbol{\textsl{k}}-1}(x)\), which follows from [18], we obtain formula ?? for \(Q_0^{\boldsymbol{\textsl{k}}}(x,y)\).
Comparing the coefficients of \(z^n\) in 33 gives the recursion identity ?? . The parity conditions for \(Q_n^{\boldsymbol{\textsl{k}}}(x,y)\)
follow from ?? and Lemma 1(iii). \(\sqcap\kern-8.0pt\sqcup\)
Proof of Proposition 2: The proof is based on the following short-time expansion of the heat kernel. This result can be found in [11], see also [12], [19]. Let \(p(t,x,y)\) be the heat kernel of the operator \(L=-\partial_x^2+V(x)\), where \(V : {\mathbb{R}}\mapsto
{\mathbb{R}}\) is a smooth function. As \(t\to 0^+\), we have the asymptotic expansion \[\label{Hadamard95expansion}
p(t,x,y)= \frac{1}{\sqrt{4\pi t}}
\exp\Big( -\frac{(y-x)^2}{4t} \Big)
\big[ 1 + u_1(x,y) t + u_2(x,y) t^2+O(t^3) \big],\tag{34}\] where \[\label{eqn95u1}
u_1(x,y)=-\int_0^1 V(x+s(y-x)) {\textrm d}s,\tag{35}\] and \[\label{eqn95u2}
u_2(x,y)=\frac{1}{2} u_1(x,y)^2 - \int_0^1
s(1-s) V''(x+s(y-x)) {\textrm d}s.\tag{36}\] The heat kernel for the operator \({\mathcal{L}}^{(1)}=-\partial_x^2+x^2\) is \[p^{(1)}(t,x,y)=\frac{ e^{-t}}{\sqrt{\pi}}
{\mathcal{E}}(e^{-2t},x,y),\] since \({\mathcal{L}}^{(1)} \psi_n=(2n+1) \psi_n\) and the functions \(\psi_n\) form an orthogonal basis for \(L_2({\mathbb{R}},{\textrm d}x)\) with \(\lVert \psi_n \rVert^2=\sqrt{\pi} 2^n n!\). Let \(\omega(x):=\partial_x^2 \ln(H_{\boldsymbol{\textsl{k}}}(x))\). From 4 we find that the heat kernel for the operator \[{\mathcal{L}}^{(2)}=-\partial_x^2+U_{\boldsymbol{\textsl{k}}}(x)+2l=-\partial_x^2+x^2-2\omega(x)+2l,\] is \[p^{(2)}(t,x,y)=\frac{e^{-t} {\mathcal{E}}_{\boldsymbol{\textsl{k}}}(e^{-2t},x,y)}{\sqrt{\pi} \psi_{\boldsymbol{\textsl{k}}}(x) \psi_{\boldsymbol{\textsl{k}}}(y)}.\] From ?? and 9 we conclude
that \[\label{p95ratio1}
\frac{p^{(2)}(t,x,y)}{p^{(1)}(t,x,y)}=\frac{1}{H_{\boldsymbol{\textsl{k}}}(x) H_{\boldsymbol{\textsl{k}}}(y)}
\sum\limits_{n=0}^{k_l+1} e^{-2nt} Q_n^{\boldsymbol{\textsl{k}}}(x,y).\tag{37}\] Let \(\alpha_1(x,y)\) and \(\alpha_2(x,y)\) (respectively, \(\beta_1(x,y)\) and \(\beta_2(x,y)\)) be the coefficients \(u_1(x,y)\) and \(u_2(x,y)\) defined by equations 35 and 36 with \(V(x)=x^2+2l-2\omega(x)\) (respectively, \(V(x)=x^2\)). For simplicity, we write \(\alpha_i=\alpha_i(x,y)\) and similarly for \(\beta_i\). According to 34 , we have \[\label{p95ratio2}
\frac{p^{(2)}(t,x,y)}{p^{(1)}(t,x,y)}=
\frac{1+\alpha_1 t+\alpha_2 t^2+O(t^3)}{1+\beta_1 t + \beta_2 t^2 + O(t^3)}.\tag{38}\]
We compute \(\beta_i\) using formulas 35 and 36 with \(V(x)=x^2\): \[\beta_1=-\frac{1}{3} (x^2+xy+y^2), \;\;\; \beta_2=\frac{1}{2} \beta_1^2 - \frac{1}{3}.\]
Next, we use the fact that \(\omega(x)=\partial_x^2 \ln(H_{\boldsymbol{\textsl{k}}}(x))\) and compute \[\Omega(x,y):=\int_0^1 \omega(x+s(y-x)) {\textrm d}s=\frac{1}{y-x} \Big( \frac{H'_{\boldsymbol{\textsl{k}}}(y)}{H_{\boldsymbol{\textsl{k}}}(y)} - \frac{H'_{\boldsymbol{\textsl{k}}}(x)}{H_{\boldsymbol{\textsl{k}}}(x)} \Big).\] After integrating by parts twice, we arrive at \[\Theta(x,y):=\int_0^1 s(1-s)\omega''(x+s(y-x)) {\textrm d}s= \frac{1}{(y-x)^2} \Big[\omega(x)+\omega(y)-2\Omega(x,y)\Big].\] Again, in what follows, we will write \(\Omega=\Omega(x,y)\) and \(\Theta=\Theta(x,y)\). Using formulas 35 and 36 with \(V(x)=x^2+2l-2\omega(x)\) we find \[\alpha_1=\beta_1-2l+2 \Omega, \;\;\; \alpha_2=\frac{1}{2} \alpha_1^2 - \frac{1}{3} + 2 \Theta.\]
We check that \[\frac{1+\alpha_1 t+\alpha_2 t^2+O(t^3)}{1+\beta_1 t + \beta_2 t^2 + O(t^3)}=
1+\gamma_1 t + \gamma_2 t^2 + O(t^3),\] where \[\gamma_1=\alpha_1-\beta_1=2 (\Omega-l),\] and \[\gamma_2=\alpha_2-\beta_2-\gamma_1 \beta_1 = 2 (\Omega-l)^2 + 2\Theta.\] Combining
the above formulas with 37 and 38 , we arrive at the following result: for all \(x,y\in{\mathbb{R}}\), as \(t\to 0^+\), \[\frac{1}{H_{\boldsymbol{\textsl{k}}}(x) H_{\boldsymbol{\textsl{k}}}(y)}
\sum\limits_{n=0}^{k_l+1} e^{-2nt} Q_n^{\boldsymbol{\textsl{k}}}(x,y)=1+2 (\Omega(x,y)-l) t + 2 \Big[ (\Omega(x,y)-l)^2 + \Theta(x,y)\Big] t^2 + O(t^3).\] Comparing the first three coefficients in the Taylor expansions at \(t=0\) on both sides of the above equation gives us the desired formulas ?? , ?? and ?? . \(\sqcap\kern-8.0pt\sqcup\)
Remark 3: Note that the asymptotic expansion 38 implies the assertion of Corollary 2, which was used crucially in the proof of
Theorem 1. Thus, when \(\boldsymbol{\textsl{k}}\) is a Krein-Adler sequence, Corollary 2 is a simple consequence of the short-time heat-kernel expansion.
Proof of Proposition 3: Let \(Q_n^k(x,y)\) be the polynomials defined in ?? . Since in this case \(\boldsymbol{\textsl{k}}=(k)\), we will write \(Q_n^k(x,y)\) instead of \(Q_n^{\boldsymbol{\textsl{k}}}(x,y)\). Similarly, we will write \(H_{k,n}(x)=\textrm{Wr}[H_k(x),H_n(x)]\) instead of \(H_{\boldsymbol{\textsl{k}},n}(x)\). We will show that \(L_n^k(x,y)=R_n^k(x,y)\) for all \(k\geqslant 0\), \(0\leqslant n \leqslant k+1\) and all \(x,y\in {\mathbb{C}}\), where \[L_n^k(x,y):=\sum_{m=0}^{n} \frac{H_m(x)H_m(y)}{2^m m!} Q_{n-m}^k(x,y), \;\;\; R_n^k(x,y):= \mathbf{1}_{\{n\ne k\}} \frac{H_{k,n}(x)H_{k,n}(y)}{2^{n+1}n!(n-k)}.\] According to Proposition 1, this implies that \(Q_n^k(x,y)\) are the polynomials appearing in ?? .
We denote \[h_m(x,y):=H_m(x)H_m(y).\] To simplify notation, we write \(h_m=h_m(x,y)\) and similarly for \(Q_n^k\), \(L_n^k\), and \(R_n^k\). From equation ?? , it follows that for \(n\geqslant 1\) and \(k\geqslant 1\), \[Q_n^k={\mathbf{1}}_{\{n\geqslant 1\}} 2k Q_{n-1}^{k-1} + {\mathbf{1}}_{\{n=1\}} h_k-2k {\mathbf{1}}_{\{n=0\}} h_{k-1}.\] Thus, for \(n\geqslant 1\) and \(k\geqslant 1\), \[L_n^k=\sum_{m=0}^{n} \bigg[ \frac{h_m}{2^m m!}\bigg] Q_{n-m}^k=2 k \sum_{m=0}^{n-1} \bigg[ \frac{h_m}{2^m m!}\bigg] Q_{n-1-m}^{k-1}+\bigg[\frac{h_{n-1}}{2^{n-1} (n-1)!}\bigg] h_k -2k \bigg[\frac{h_{n}}{2^{n} n!}\bigg] h_{k-1},\] and we obtain the recurrence relation \[\label{L95recursion} L_n^k=2 k L_{n-1}^{k-1}+\frac{nh_k h_{n-1}-k h_{k-1} h_n}{2^{n-1} n!}.\tag{39}\]
Next, using the facts \[H_n'(x)=2n H_{n-1}(x)=2x H_n(x)-H_{n+1}(x),\] we check that \[\begin{align} &H_{k,n}(x)=2n H_k(x) H_{n-1}(x)-2k H_{k-1}(x)H_n(x), \\ &H_{k-1,n-1}(x)=H_{k}(x)H_{n-1}(x) -H_{k-1}(x) H_{n}(x). \end{align}\] Using the above two formulas and a straightforward but tedious calculation, we verify that \[H_{k,n}(x)H_{k,n}(y)=4n(n-k)h_k h_{n-1} - 4 k (n-k) h_{k-1} h_n+4 n k H_{k-1,n-1}(x) H_{k-1,n-1}(y).\] When \(n\geqslant 1\) and \(n\neq k\), we divide both sides of the above identity by \(2^{n+1} n! (n-k)\) and obtain \[\label{R95recursion} R_n^k=2 k R_{n-1}^{k-1}+\frac{nh_k h_{n-1}-k h_{k-1} h_n}{2^{n-1} n!}.\tag{40}\] We check that the above identity also holds for \(n=k\), since in this case \(R_n^k=R_{n-1}^{k-1}=nh_k h_{n-1}-k h_{k-1} h_n=0\).
We see from formulas 39 and 40 that \(L_n^k\) and \(R_n^k\) satisfy the same recurrence relation, which lowers both indices
\(n\) and \(k\) by one. We want to prove the identity \(L_n^k=R_n^k\) for \(k\geqslant 0\) and \(0\leqslant n \leqslant k+1\). After applying the recurrence formulas 39 and 40 \(M=\min(n,k)\) times, we arrive at either the case
\(n=0\), \(k\geqslant 0\) or \(n=1\), \(k=0\). In other words, if the identity \(L_n^k=R_n^k\) holds for \(n=0\), \(k\geqslant 0\) and \(n=1\), \(k=0\), then it holds for all
\(k\geqslant 0\) and \(n=0,1,\dots,k+1\). The fact that the identity \(L_n^k=R_n^k\) holds in these boundary cases is easily verified directly using the
formulas \[H_{k,0}(x)=-2k H_{k-1}(x), \;\;\; H_{0,n}(x)=2n H_{n-1}(x),\] and the definition ?? of \(Q_n^k(x,y)\). \(\sqcap\kern-8.0pt\sqcup\)
Dept. of Mathematics and Statistics, York University, 4700 Keele Street, Toronto, ON, M3J 1P3, Canada.
Email: akuznets@yorku.ca, yuanm@yorku.ca↩︎