Multi-variable special polynomials using an operator ordering method

Abstract

Using an operator ordering method for some commutative superposition operators, we introduce two new multi-variable special polynomials and their generating functions, and present some new operator identities and integral formulas involving the two special polynomials. Instead of calculating complicated partial differential, we use the special polynomials and their generating functions to concisely address the normalization, photocount distributions and Wigner distributions of several quantum states that can be realized physically, the results of which provide real convenience for further investigating the properties and applications of these states.

Graphical abstract

Keywords

multi-variable special polynomial / generating function / operator ordering method / new operator identity and integral formula / Wigner function

1 Introduction

As a class of typical special polynomials, Hermite polynomials play a prominent role in various fields of mathematics and physics due to their interesting basic properties (e.g., orthogonality and completeness) and some important relevant identities such as the generating function, product of polynomials, recurrence relation, and differential equation.

The one-variable Hermite polynomials H_n(x) are usually defined via their generating functions [1, 2]

\begin{array}{}\text{(1)}& {\text{e}}^{-t^2+2xt}={\displaystyle \sum _{n=0}^{\infty }}\frac{t^n}{n!}{H}_n(x)\mathrm{.}\end{array}

H_n(x) are eigenstates of the quantum Harmonic oscillator and eigenfunctions of the fractional Fourier transformation [3]. Also, they can help to solve the operator Fredholm equation [4] and the eigenvalue problems of coupled harmonic oscillator systems and generalized angular momentum systems [5, 6]. In the two-variable case, Hermite polynomials H_{n\mathrm{,}m}(x\mathrm{,}y) can be defined via [7-9]

\begin{array}{}\text{(2)}& {\text{e}}^{-t{t}^{\prime }+tx+t^{\prime }y}={\displaystyle \sum _{m\mathrm{,}n=0}^{\infty }}\frac{t^n{t^{\prime }}^m}{n!m!}{H}_{n\mathrm{,}m}(x\mathrm{,}y)\mathrm{,}\end{array}

and their partial differential representations and power-series expansions are respectively obtained as

\begin{array}{}\text{(3)}& \begin{array}{ll}{H}_{n\mathrm{,}m}(x\mathrm{,}y)\hfill & =\frac{{\partial }^{n+m}}{\partial t^n\partial t^{\text{'}m}}{\text{e}}^{-t{t}^{\prime }+tx+t^{\prime }y}{|}_{t=t^{\prime }=0}\hfill \\ \hfill & ={\displaystyle \sum _{l=0}^{\min (m\mathrm{,}n)}}\left(\begin{array}{c}n\\ l\end{array}\right)\left(\begin{array}{c}m\\ l\end{array}\right)l!{\left(-1\right)}^l{x}^{n-l}{y}^{m-l}\mathrm{.}\hfill \end{array}\end{array}

H_{n\mathrm{,}m}(x\mathrm{,}y) can be explained as the transition amplitudes of number states in the dynamics of the forced Harmonic oscillator [10] or the eigenfunctions of the two-dimensional complex fractional Fourier transform [11], and they are useful in studying the Talbot effect in a quadratic-index medium [12], the Bargmann transformation, and the quantum entanglement phenomena. Besides, both the one- and two-variable Hermite polynomials can be applied to deduce some new asymptotic formulas [13] and bosonic operator ordering identities [14], and produce some new non-Gaussian quantum states [15-18] that may be used as a key entangled resource for realizing certain quantum information tasks including quantum teleportation, metrology and communications. Hence, up to now, many kinds of deformed Hermite polynomials, such as degenerate Hermite polynomials [2], holomorphic Hermite polynomials [19], and q-Hermite polynomials [20], are successively proposed and widely used in the probability theory, graph theory, number theory, and other areas of mathematics and physics.

Motivated by the universality and importance for applications of Hermite polynomials in the physical realm, in this article we propose two new multi-variable special polynomials via generalizing the power-series expansions of the ordinary Hermite polynomials H_{n\mathrm{,}m}(x\mathrm{,}y). The generalization mainly concerns the terms {(-1)}^l and x^{n-l}{y}^{m-l} in Eq. (3), which would be respectively replaced by a more general power function {\vartheta }^l (\vartheta being an arbitrary parameter) and the product of Hermite polynomials H_{n-l}\left(x/2\right)H_{m-l}\left(y/2\right), leading to the generalized form

\begin{array}{}\text{(4)}& {\displaystyle \sum _{l=0}^{\min (n\mathrm{,}m)}}\left(\begin{array}{c}n\\ l\end{array}\right)\left(\begin{array}{c}m\\ l\end{array}\right)l!{\vartheta }^l{H}_{n-l}\left(\frac{x}{2}\right)H_{m-l}\left(\frac{y}{2}\right)\mathrm{.}\end{array}

Given Eq. (4), two interesting questions naturally arise: can a new useful special polynomial appear? If yes, what are its applications in quantum optics? Furthermore, what would happen if the term x^{n-l}{y}^{m-l} is replaced by the product of two-variable Hermite polynomials H_{n-i\mathrm{,}m-j}\left(x\mathrm{,}y\right)H_{l-i\mathrm{,}k-j}\left(x^{\prime }\mathrm{,}{y}^{\prime }\right) rather than H_{n-l}\left(x/2\right)H_{m-l}\left(y/2\right)? To tackle these problems and avoid the complications arising from the summations of the products of Hermite polynomials, we shall fully use an operator ordering method for some commutative superposition operators. Indeed, the operator ordering method has made a significant contribution to the development of quantum physics in many aspects, such as the operator integration theory, operator special polynomial theory, and representation theory [21, 22].

The rest of the article is organized as follows. In Section 2, we introduce two new multi-variable special polynomials and their generating functions using the operator ordering method. In Section 3, we derive some new operator identities and integral formulas of the two new special polynomials. Applications of the two special polynomials in the normalization, photocount distributions and Wigner distributions of several quantum states that can be physically implemented are discussed in Section 4. The main results of this article are summarized in Section 5.

2 Multi-variable special polynomials and their generating functions

We first examine how to use the operator ordering method for some commutative superposition operators to introduce two new multi-variable special polynomials and their generating functions.

2.1 Three-variable case

In terms of the definition of H_{n\mathrm{,}m}(x\mathrm{,}y), we expand the polynomials {\mathcal{H}}_{n\mathrm{,}m}(x\mathrm{,}y;\vartheta ) as

\begin{array}{}\text{(5)}& \begin{array}{l}{\mathcal{H}}_{n\mathrm{,}m}(x\mathrm{,}y;\vartheta )\\ =\frac{{\partial }^{n+m}}{\partial s^n\partial {\tau }^m}{\exp \left(\vartheta s\tau +sx+\tau y\right)|}_{s=\tau =0}\\ ={\displaystyle \sum _{l=0}^{\min (n\mathrm{,}m)}}\left(\begin{array}{c}n\\ l\end{array}\right)\left(\begin{array}{c}m\\ l\end{array}\right)l!{\vartheta }^l{x}^{n-l}{y}^{m-l}\mathrm{,}\end{array}\end{array}

which are indeed two-variable Hermite polynomials since the polynomials {\mathcal{H}}_{n\mathrm{,}m}(x\mathrm{,}y;\vartheta )reduce to H_{n\mathrm{,}m}(x\mathrm{,}y) via replacing i\sqrt{\vartheta }s\to s^{\text{'}} and i\sqrt{\vartheta }\tau \to {\tau }^{\text{'}} in Eq. (5). Supposing that X=\sqrt{2}\left(a+a^{†}\right) and Y=\sqrt{2}\left(b+b^{†}\right), where a^{†}\mathrm{,}{b}^{†} are respectively the creation operators of two boson modes, and noting that \left[X\mathrm{,}Y\right]=0, thus we make the substitutions x\to Xand y\to Y in Eq. (5), and introduce the operator identity

\begin{array}{}\text{(6)}& \begin{array}{l}{\mathcal{H}}_{n\mathrm{,}m}(X\mathrm{,}Y;\vartheta )\\ =\frac{{\partial }^{n+m}}{\partial s^n\partial {\tau }^m}{\exp \left(\vartheta s\tau +sX+\tau Y\right)|}_{s=\tau =0}\\ ={\displaystyle \sum _{l=0}^{\min (n\mathrm{,}m)}}\left(\begin{array}{c}n\\ l\end{array}\right)\left(\begin{array}{c}m\\ l\end{array}\right)l!{\vartheta }^l{X}^{n-l}{Y}^{m-l}\mathrm{.}\end{array}\end{array}

Using the following anti-normal ordering products (indicated by the symbol ⋮⋮ [23-25])

\begin{array}{}\text{(7)}& X^n=⋮H_n\left(\frac{X}{2}\right)⋮\mathrm{,}{Y}^m=⋮H_m\left(\frac{Y}{2}\right)⋮\mathrm{,}\end{array}

we rewrite Eq. (6) as

\begin{array}{}\text{(8)}& \begin{array}{ll}{\mathcal{H}}_{n\mathrm{,}m}(X\mathrm{,}Y;\vartheta )\hfill & ={\displaystyle \sum _{l=0}^{\min (n\mathrm{,}m)}}\left(\begin{array}{c}n\\ l\end{array}\right)\left(\begin{array}{c}m\\ l\end{array}\right)l!{\vartheta }^l\hfill \\ \hfill & \times ⋮H_{n-l}\left(\frac{X}{2}\right)H_{m-l}\left(\frac{Y}{2}\right)⋮\mathrm{.}\hfill \end{array}\end{array}

On the other hand, using the Baker–Hausdorff formula, i,e., {\text{e}}^{A+B}={\text{e}}^A{\text{e}}^B{\text{e}}^{-[A\mathrm{,}B]/2}={\text{e}}^B{\text{e}}^A{\text{e}}^{[A\mathrm{,}B]/2}, which holds for \left[\left[A\mathrm{,}B\right]\mathrm{,}A\right]=\left[\left[A\mathrm{,}B\right]\mathrm{,}B\right]=0, we give the operator identity as

\begin{array}{}\text{(9)}& \begin{array}{l}\exp \left(\vartheta s\tau +sX+\tau Y\right)\\ =⋮\exp \left(sX+\tau Y+\vartheta s\tau -s^2-{\tau }^2\right)⋮\mathrm{.}\end{array}\end{array}

Thus, combining Eqs. (6), (8) and (9) leads to

\begin{array}{}\text{(10)}& \begin{array}{l}\frac{{\partial }^{n+m}}{\partial s^n\partial {\tau }^m}{⋮\exp \left(sX+\tau Y+\vartheta s\tau -s^2-{\tau }^2\right)⋮|}_{s=\tau =0}\\ ={\displaystyle \sum _{l=0}^{\min (n\mathrm{,}m)}}\left(\begin{array}{c}n\\ l\end{array}\right)\left(\begin{array}{c}m\\ l\end{array}\right)l!{\vartheta }^l⋮H_{n-l}\left(\frac{X}{2}\right)H_{m-l}\left(\frac{Y}{2}\right)⋮\mathrm{.}\end{array}\end{array}

Observing that both sides of Eq. (10) are in anti-normal ordering, thus replacing X\to xand Y\to y in Eq. (10) and comparing with Eq. (3), we naturally obtain

\begin{array}{}\text{(11)}& \begin{array}{l}\exp \left(-s^2-{\tau }^2+\vartheta s\tau +sx+\tau y\right)\\ ={\displaystyle \sum _{n\mathrm{,}m=0}^{\infty }}\frac{s^n{\tau }^m}{n!m!}{ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right)\mathrm{,}\end{array}\end{array}

where

\begin{array}{}\text{(12)}& \begin{array}{l}{ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right)\\ =\frac{{\partial }^{n+m}}{\partial s^n\partial {\tau }^m}{\exp \left(sx+\tau y+\vartheta s\tau -s^2-{\tau }^2\right)|}_{s=\tau =0}\\ ={\displaystyle \sum _{l=0}^{\min (n\mathrm{,}m)}}\left(\begin{array}{c}n\\ l\end{array}\right)\left(\begin{array}{c}m\\ l\end{array}\right)l!{\vartheta }^l{H}_{n-l}\left(\frac{x}{2}\right)H_{m-l}\left(\frac{y}{2}\right)\end{array}\end{array}

is a new two-index, three-variable special polynomial with generating function {\text{e}}^{-s^2-{\tau }^2+\vartheta s\tau +sx+\tau y}. Specially, when \vartheta =0, {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right) becomes the product of Hermite polynomials H_n\left(x/2\right)H_m\left(y/2\right). Further, using the differential identity H_m^{\prime }\left(x\right)=2m{H}_{m-1}\left(x\right), we arrive at the first-order partial differential equations of {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right) with respect to the variables x and y, i.e.,

\begin{array}{}\text{(13)}& \begin{array}{l}\frac{\partial }{\partial x}{ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right)=n{ℌ}_{n-\mathrm{1,}m}\left(x\mathrm{,}y;\vartheta \right)\mathrm{,}\\ \frac{\partial }{\partial y}{ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right)=m{ℌ}_{n\mathrm{,}m-1}\left(x\mathrm{,}y;\vartheta \right)\mathrm{.}\end{array}\end{array}

Thus, its high-order differential identity yields

\begin{array}{}\text{(14)}& \begin{array}{l}\frac{{\partial }^{k+l}}{\partial x^k\partial y^l}{ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right)\\ =\frac{n!m!}{\left(n-k\right)!\left(m-l\right)!}{ℌ}_{n-k\mathrm{,}m-l}\left(x\mathrm{,}y;\vartheta \right)\mathrm{,}\end{array}\end{array}

which has the same form as the well-known differential of the ordinary Hermite polynomials H_{n\mathrm{,}m}\left(x\mathrm{,}y\right).

2.2 Six-variable case

To obtain six-variable special polynomials, we first introduce the product of Hermite polynomials {\mathcal{H}}_{n\mathrm{,}l}(x\mathrm{,}{x}^{\prime };\nu ){\mathcal{H}}_{m\mathrm{,}k}(y\mathrm{,}{y}^{\prime };\upsilon ) and denote it as {\mathcal{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}\left(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon \right) that is expressed as

\begin{array}{}\text{(15)}& \begin{array}{l}{\mathcal{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}\left(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon \right)\\ =\frac{{\partial }^{m+n}}{\partial s^n\partial {\tau }^m}\frac{{\partial }^{l+k}}{\partial {s^{\prime }}^l\partial {{\tau }^{\prime }}^k}\exp (\nu s{s}^{\prime }+\upsilon \tau {\tau }^{\prime }\\ {+sx+s^{\prime }{x}^{\prime }+\tau y+{\tau }^{\prime }{y}^{\prime })|}_{s=s^{\prime }=\tau ={\tau }^{\prime }=0}\mathrm{,}\end{array}\end{array}

where \nu \mathrm{,}\upsilon are taken as arbitrary parameters. Introducing the superposed operators W=a+b^{†}, Z=a^{†}+b, W^{\prime }=c+d^{†} and Z^{\prime }=c^{†}+d, where a^{†}\mathrm{,}{b}^{†}\mathrm{,}{c}^{†}\mathrm{,}{d}^{†} are respectively the creation operators of four bosonic modes, and noting that the operators W\mathrm{,}Z\mathrm{,}{W}^{\prime } and Z^{\prime } commute with each other, so we make the substitutions x\to W, y\to Z, x^{\prime }\to W^{\prime } and y^{\prime }\to Z^{\prime } in Eq. (15) and obtain the operator identity as

\begin{array}{}\text{(16)}& \begin{array}{l}{\mathcal{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(W\mathrm{,}Z\mathrm{,}{W}^{\prime }\mathrm{,}{Z}^{\prime };\nu \mathrm{,}\upsilon )\\ =\frac{{\partial }^{m+n}}{\partial s^n\partial {\tau }^m}\frac{{\partial }^{l+k}}{\partial {s^{\prime }}^l\partial {{\tau }^{\prime }}^k}\exp (\nu s{s}^{\prime }+\upsilon \tau {\tau }^{\prime }\\ {+sW+s^{\prime }{W}^{\prime }+\tau Z+{\tau }^{\prime }{Z}^{\prime })|}_{s=s^{\prime }=\tau ={\tau }^{\prime }=0}\\ ={\displaystyle \sum _{i\mathrm{.}j=0}^{\min (n\mathrm{,}m\mathrm{,}l\mathrm{,}k)}}\left(\begin{array}{c}n\\ i\end{array}\right)\left(\begin{array}{c}m\\ j\end{array}\right)\left(\begin{array}{c}l\\ i\end{array}\right)\left(\begin{array}{c}k\\ j\end{array}\right)\\ \times i!j!{\nu }^i{\upsilon }^j{W}^{n-i}{W^{\prime }}^{l-i}{Z}^{m-j}{Z^{\prime }}^{k-j}\mathrm{.}\end{array}\end{array}

Further, using the operator identities

\begin{array}{}\text{(17)}& W^n{Z}^m=⋮H_{n\mathrm{,}m}\left(W\mathrm{,}Z\right)⋮\mathrm{,}{W^{\prime }}^l{Z^{\prime }}^k=⋮H_{l\mathrm{,}k}\left(W^{\prime }\mathrm{,}{Z}^{\prime }\right)⋮\mathrm{,}\end{array}

which are given by the s-ordered expansions of operators and the completeness relation of the entangled states that are the eigenstates of the operator pair (W\mathrm{,}Z) or (W^{\prime }\mathrm{,}{Z}^{\prime }) [22, 26], thus Eq. (16) becomes

\begin{array}{}\text{(18)}& \begin{array}{l}{\mathcal{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(W\mathrm{,}Z\mathrm{,}{W}^{\prime }\mathrm{,}{Z}^{\prime };\nu \mathrm{,}\upsilon )\\ ={\displaystyle \sum _{i\mathrm{.}j=0}^{\min (n\mathrm{,}m\mathrm{,}l\mathrm{,}k)}}\left(\begin{array}{c}n\\ i\end{array}\right)\left(\begin{array}{c}m\\ j\end{array}\right)\left(\begin{array}{c}l\\ i\end{array}\right)\left(\begin{array}{c}k\\ j\end{array}\right)i!j!\\ \times {\nu }^i{\upsilon }^j⋮H_{n-i\mathrm{,}m-j}\left(W\mathrm{,}Z\right)H_{l-i\mathrm{,}k-j}\left(W^{\prime }\mathrm{,}{Z}^{\prime }\right)⋮\mathrm{,}\end{array}\end{array}

which is just the antinormal ordering product of {\mathcal{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(W\mathrm{,}Z\mathrm{,}{W}^{\prime }\mathrm{,}{Z}^{\prime };\nu \mathrm{,}\upsilon ). Besides, using the Baker–Hausdorff formula, we have

\begin{array}{}\text{(19)}& \begin{array}{l}\exp \left(\nu s{s}^{\prime }+\upsilon \tau {\tau }^{\prime }+sW+s^{\prime }{W}^{\prime }+\tau Z+{\tau }^{\prime }{Z}^{\prime }\right)\\ =⋮\exp (-s\tau -s^{\prime }{\tau }^{\prime }+\nu s{s}^{\prime }+\upsilon \tau {\tau }^{\prime }+sW\\ +s^{\prime }{W}^{\prime }+\tau Z+{\tau }^{\prime }{Z}^{\prime })⋮\mathrm{.}\end{array}\end{array}

Comparing with Eqs. (16), (18) and (19), and replacing W\to x\mathrm{,}Z\to y\mathrm{,}{W}^{\prime }\to x^{\prime } and Z^{\prime }\to y^{\prime } in the result after comparison, we obtain a new four-index, six-variable special polynomial as

\begin{array}{}\text{(20)}& \begin{array}{l}\frac{{\partial }^{m+n}}{\partial s^n\partial {\tau }^m}\frac{{\partial }^{l+k}}{\partial {s^{\prime }}^l\partial {{\tau }^{\prime }}^k}\exp (-s\tau -s^{\prime }{\tau }^{\prime }+\nu s{s}^{\prime }\\ {+\upsilon \tau {\tau }^{\prime }+xs+\sout{x}{s}^{\prime }+y\tau +y^{\prime }{\tau }^{\prime })|}_{s=s^{\prime }=\tau ={\tau }^{\prime }=0}\\ ={\displaystyle \sum _{i\mathrm{.}j=0}^{\min (n\mathrm{,}m\mathrm{,}l\mathrm{,}k)}}\left(\begin{array}{c}n\\ i\end{array}\right)\left(\begin{array}{c}m\\ j\end{array}\right)\left(\begin{array}{c}l\\ i\end{array}\right)\left(\begin{array}{c}k\\ j\end{array}\right)i!j!\\ \times {\nu }^i{\upsilon }^j{H}_{n-i\mathrm{,}m-j}\left(x\mathrm{,}y\right)H_{l-i\mathrm{,}k-j}\left(x^{\prime }\mathrm{,}{y}^{\prime }\right)\\ \equiv {\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon )\mathrm{,}\end{array}\end{array}

such that its generating function is

\begin{array}{}\text{(21)}& \begin{array}{l}\exp (-s\tau -s^{\prime }{\tau }^{\prime }+\nu s{s}^{\prime }+\upsilon \tau {\tau }^{\prime }+xs+x^{\prime }{s}^{\prime }+y\tau +y^{\prime }{\tau }^{\prime })\\ ={\displaystyle \sum _{n\mathrm{,}m\mathrm{,}l\mathrm{,}k=0}^{\infty }}\frac{s^n{\tau }^m{s^{\prime }}^l{{\tau }^{\prime }}^k}{n!m!l!k!}{\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon )\mathrm{.}\end{array}\end{array}

Similarly, using the differential identities of Hermite polynomials H_{n\mathrm{,}m}\left(x\mathrm{,}y\right), we also arrive at the high-order partial differential of {\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon ) with respect to the variables x\mathrm{,}y\mathrm{,}{x}^{\prime } and y^{\prime }, i.e.,

\begin{array}{}\text{(22)}& \begin{array}{l}\frac{{\partial }^{i+j}}{\partial x^i\partial y^j}\frac{{\partial }^{i^{\prime }+j^{\prime }}}{\partial {x^{\prime }}^{i^{\prime }}\partial {y^{\prime }}^{j^{\prime }}}{\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon )\\ =\frac{n!m!l!k!}{\left(n-i\right)!\left(m-j\right)!(l-i^{\prime })!(k-j^{\prime })!}\\ \times {\mathfrak{F}}_{n-i\mathrm{,}m-j\mathrm{,}l-i^{\text{'}}\mathrm{,}k-j^{\text{'}}}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon )\mathrm{.}\end{array}\end{array}

3 New operator identities and integral formulas

In this section, we show how the two new special polynomials and their generating functions can be used to deduce some new operator identities and integral formulas.

By comparing of Eqs. (8) and (12), we obtain a new operator identity

\begin{array}{}\text{(23)}& {\mathcal{H}}_{n\mathrm{,}m}(X\mathrm{,}Y;\vartheta )=⋮{ℌ}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right)⋮\mathrm{,}\end{array}

which is indeed useful to prove many interesting relations about {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right). For instance, from the generating function of {\mathcal{H}}_{n\mathrm{,}m}(X\mathrm{,}Y;\vartheta ), we obtain the recurrence relation of {\mathcal{H}}_{n\mathrm{,}m}(X\mathrm{,}Y;\vartheta ) as

\begin{array}{}\text{(24)}& \begin{array}{c}nm{\mathcal{H}}_{n-\mathrm{1,}m-1}\left(X\mathrm{,}Y;\vartheta \right)-nX{\mathcal{H}}_{n-\mathrm{1,}m}\left(X\mathrm{,}Y;\vartheta \right)\\ +n{\mathcal{H}}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right)=0.\end{array}\end{array}

Substituting Eq. (23) into Eq. (24) leads to the form

\begin{array}{}\text{(25)}& \begin{array}{c}nm⋮{ℌ}_{n-\mathrm{1,}m-1}\left(X\mathrm{,}Y;\vartheta \right)⋮-nX⋮{ℌ}_{n-\mathrm{1,}m}\left(X\mathrm{,}Y;\vartheta \right)⋮\\ +n⋮{ℌ}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right)⋮=\mathrm{0,}\end{array}\end{array}

such that the recurrence relation of {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right) is

\begin{array}{}\text{(26)}& \begin{array}{c}nm{ℌ}_{n-\mathrm{1,}m-1}\left(x\mathrm{,}y;\vartheta \right)-nx{ℌ}_{n-\mathrm{1,}m}\left(x\mathrm{,}y;\vartheta \right)\\ +n{ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right)=\mathrm{0,}\end{array}\end{array}

which is the same as that of the ordinary Hermite polynomials H_{n\mathrm{,}m}\left(x\mathrm{,}y\right) in form.

From Eqs. (6) and (11), we have

\begin{array}{}\text{(27)}& \begin{array}{l}{\displaystyle \sum _{n\mathrm{,}m=0}^{\infty }}\frac{s^n{\tau }^m}{n!m!}{ℌ}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right)\\ ={\text{e}}^{-s^2-{\tau }^2+\vartheta s\tau +sX+\tau Y}=:{\text{e}}^{\vartheta s\tau +sX+\tau Y}:\\ ={\displaystyle \sum _{n\mathrm{,}m=0}^{\infty }}\frac{s^n{\tau }^m}{n!m!}:{\mathcal{H}}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right):\mathrm{,}\end{array}\end{array}

from which it then follows that

\begin{array}{}\text{(28)}& {ℌ}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right)=:{\mathcal{H}}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right):\mathrm{,}\end{array}

where :: denotes normal ordering [27-29]. Operating both sides of Eq. (28) on two-mode vacuum, we obtain a new quantum state as {\mathcal{H}}_{n\mathrm{,}m}\left(\sqrt{2}{a}^{†}\mathrm{,}\sqrt{2}{b}^{†};\vartheta \right)|00\rangle. Using Eq. (28) and the completeness of coordinate representation in normal ordering, we obtain

\begin{array}{}\text{(29)}& \begin{array}{l}{ℌ}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right)\\ =\frac{1}{\pi }{\displaystyle {\int }_{-\infty }^{\infty }}\text{d}{q}_1\text{d}{q}_2:{\text{e}}^{-{(q_1-Q_1)}^2-{(q_2-Q_2)}^2}:{ℌ}_{n\mathrm{,}m}(2q_1\mathrm{,2}{q}_2;\vartheta )\\ =:{\mathcal{H}}_{n\mathrm{,}m}\left(X\mathrm{,}Y;\vartheta \right):\mathrm{,}\end{array}\end{array}

from which we have a new integral formula

\begin{array}{}\text{(30)}& \begin{array}{l}\frac{1}{\pi }{\displaystyle {\int }_{-\infty }^{\infty }}\text{d}{q}_1\text{d}{q}_2{\text{e}}^{-{(q_1-x)}^2-{(q_2-y)}^2}{ℌ}_{n\mathrm{,}m}(2q_1\mathrm{,2}{q}_2;\vartheta )\\ ={\mathcal{H}}_{n\mathrm{,}m}\left(2x\mathrm{,2}y;\vartheta \right)\mathrm{.}\end{array}\end{array}

On the other hand, replacing x\to \sqrt{\vartheta }W and y\to \sqrt{{\vartheta }^{\ast }}Z in Eq. (11), we have

\begin{array}{}\text{(31)}& \begin{array}{l}{\text{e}}^{-s^2-{\tau }^2+\left|\vartheta \right|s\tau +\sqrt{\vartheta }sW+\sqrt{{\vartheta }^{\ast }}\tau Z}\\ ={\displaystyle \sum _{n\mathrm{,}m=0}^{\infty }}\frac{s^n{\tau }^m}{n!m!}{ℌ}_{n\mathrm{,}m}\left(\sqrt{\vartheta }W\mathrm{,}\sqrt{{\vartheta }^{\ast }}Z;\vartheta \right)\mathrm{.}\end{array}\end{array}

Further, using the Baker-Hausdorff formula and Eq. (1), Eq. (31) becomes

\begin{array}{}\text{(32)}& \begin{array}{l}⋮{\text{e}}^{-s^2-{\tau }^2+\sqrt{\vartheta }sW+\sqrt{{\vartheta }^{\ast }}\tau Z}⋮\\ ={\displaystyle \sum _{n\mathrm{,}m=0}^{\infty }}\frac{s^n{\tau }^m}{n!m!}⋮H_n\left(\frac{\sqrt{\vartheta }W}{2}\right)H_m\left(\frac{\sqrt{{\vartheta }^{\ast }}Z}{2}\right)⋮\mathrm{.}\end{array}\end{array}

Thus, comparing between Eq. (31) and Eq. (32) leads to the operator identity

\begin{array}{}\text{(33)}& \begin{array}{l}{ℌ}_{n\mathrm{,}m}\left(\sqrt{\vartheta }W\mathrm{,}\sqrt{{\vartheta }^{\ast }}Z;\left|\vartheta \right|\right)\\ =⋮H_n\left(\frac{\sqrt{\vartheta }W}{2}\right)H_m\left(\frac{\sqrt{{\vartheta }^{\ast }}Z}{2}\right)⋮\mathrm{.}\end{array}\end{array}

Besides, according to Eqs. (1) and (11), we obtain

\begin{array}{}\text{(34)}& \begin{array}{l}{\text{e}}^{-s^2-{\tau }^2+\sqrt{\vartheta }sW+\sqrt{{\vartheta }^{\ast }}\tau Z}\\ ={\displaystyle \sum _{n\mathrm{,}m=0}^{\infty }}\frac{s^n{\tau }^m}{n!m!}{H}_n\left(\frac{\sqrt{\vartheta }W}{2}\right)H_m\left(\frac{\sqrt{{\vartheta }^{\ast }}Z}{2}\right)\\ =:{\text{e}}^{-s^2-{\tau }^2+\left|\vartheta \right|s\tau +\sqrt{\vartheta }sW+\sqrt{{\vartheta }^{\ast }}\tau Z}:\\ ={\displaystyle \sum _{n\mathrm{,}m=0}^{\infty }}\frac{s^n{\tau }^m}{n!m!}:{ℌ}_{n\mathrm{,}m}\left(\sqrt{\vartheta }W\mathrm{,}\sqrt{{\vartheta }^{\ast }}Z;\left|\vartheta \right|\right):\mathrm{,}\end{array}\end{array}

which leads to another operator identity

\begin{array}{}\text{(35)}& \begin{array}{l}{H}_n\left(\frac{\sqrt{\vartheta }W}{2}\right)H_m\left(\frac{\sqrt{{\vartheta }^{\ast }}Z}{2}\right)\\ =:{ℌ}_{n\mathrm{,}m}\left(\sqrt{\vartheta }W\mathrm{,}\sqrt{{\vartheta }^{\ast }}Z;\left|\vartheta \right|\right):\end{array}\end{array}

or another new quantum state {ℌ}_{n\mathrm{,}m}(\sqrt{\vartheta }{b}^{†}\mathrm{,}\sqrt{{\vartheta }^{\ast }}{a}^{†}; \left|\vartheta \right|)|00\rangle. Further, using completeness relation of the entangled states |\varsigma \rangle, where the state |\varsigma \rangle yields the eigenvector equations, i.e., W|\varsigma \rangle =\varsigma |\varsigma \rangle and Z|\varsigma \rangle ={\varsigma }^{\ast }|\varsigma \rangle, we have

\begin{array}{}\text{(36)}& \begin{array}{l}{H}_n\left(\frac{\sqrt{\vartheta }W}{2}\right)H_m\left(\frac{\sqrt{{\vartheta }^{\ast }}Z}{2}\right)\\ ={\displaystyle \int }\frac{{\text{d}}^2\varsigma }{\pi }:{\text{e}}^{-\left(\varsigma -W\right)\left({\varsigma }^{\ast }-Z\right)}:H_n\left(\frac{\sqrt{\vartheta }\varsigma }{2}\right)H_m\left(\frac{\sqrt{{\vartheta }^{\ast }}{\varsigma }^{\ast }}{2}\right)\\ =:{ℌ}_{n\mathrm{,}m}\left(\sqrt{\vartheta }W\mathrm{,}\sqrt{{\vartheta }^{\ast }}Z;\left|\vartheta \right|\right):\mathrm{,}\end{array}\end{array}

from which we get another new integral formula

\begin{array}{}\text{(37)}& \begin{array}{l}{\displaystyle \int }\frac{{\text{d}}^2\varsigma }{\pi }{\text{e}}^{-\left(\varsigma -\sigma \right)\left({\varsigma }^{\ast }-{\sigma }^{\ast }\right)}{H}_n\left(\frac{\sqrt{\vartheta }\varsigma }{2}\right)H_m\left(\frac{\sqrt{{\vartheta }^{\ast }}{\varsigma }^{\ast }}{2}\right)\\ ={ℌ}_{n\mathrm{,}m}\left(\sqrt{\vartheta }\sigma \mathrm{,}\sqrt{{\vartheta }^{\ast }}{\sigma }^{\ast };\left|\vartheta \right|\right)\mathrm{.}\end{array}\end{array}

Indeed, Eq. (37) has some practical uses in quantum theory, e.g. calculating the normalization of the photon-modulated vacuum state {\left(ta+r{a}^{†}\right)}^m|0\rangle, as shown below. Similarly, using the operator ordering method, we can also derive some new relevant identities involving the polynomials {\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon ). For instance, comparison of Eqs. (18–20) leads to the operator identity

\begin{array}{}\text{(38)}& \begin{array}{l}{\mathcal{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(W\mathrm{,}Z\mathrm{,}{W}^{\prime }\mathrm{,}{Z}^{\prime };\nu \mathrm{,}\upsilon )\\ =⋮{\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(W\mathrm{,}Z\mathrm{,}{W}^{\prime }\mathrm{,}{Z}^{\prime };\nu \mathrm{,}\upsilon )⋮\mathrm{.}\end{array}\end{array}

However, we will not explore them here owing to the similarity of some identities involving the two special polynomials and the limited space of this article.

4 Applications of multi-variable special polynomials

In this section, we present how to use the two new multi-variable special polynomials to more conveniently investigate some fundamental topics in quantum optics, where, in particular, the normalization, photocount distribution and Wigner distribution are emphasized.

4.1 Normalization

Normalization is significant for characterizing the success probability in a state preparation scheme and further investigating the properties and applications of this state. Here, we first use the two new special polynomials to calculate the normalization of several physically achievable quantum states. Theoretically, the state a^{m}S\left(r\right)|0\rangle can be obtained via repeatedly acting the annihilation operator a (i.e., single-photon substraction) on a squeezed vacuum S\left(r\right)|0\rangle (S\left(r\right) being single-mode squeezing operator [30, 31]) for m times [32], and the experimental implementation should be feasible for small m, since the single-photon subtraction has been achieved via the high transmittance beamsplitter [33]. Indeed, the state aS\left(r\right)|0\rangle has been successfully prepared with a periodically-poled KTiOPO₄ crystal [34]. Using the operator identity a^{†m}{a}^m=⋮H_{m\mathrm{,}m}(a^{†}\mathrm{,}a)⋮ [35], we give the normalized factor of the state a^{m}S\left(r\right)|0\rangle as

\begin{array}{}\text{(39)}& \begin{array}{ll}{N}_m\hfill & =\text{sech}r\langle 0|⋮\exp \left(\frac{1}{2}{a}^2\tanh r\right)\hfill \\ \hfill & \times H_{m\mathrm{,}m}(a^{†}\mathrm{,}a)\exp \left(\frac{1}{2}{a}^{†2}\tanh r\right)⋮|0\rangle \mathrm{.}\hfill \end{array}\end{array}

Inserting the completeness of coherent states [36] into Eq. (39), and using the generating function of H_{n\mathrm{,}m}(x\mathrm{,}y) and the integral formula

we have

\begin{array}{}\text{(41)}& N_m=\frac{{\partial }^{2m}}{\partial s^m\partial {\tau }^m}{{\text{e}}^{\frac{1}{4}({\tau }^2+s^2)\sinh 2r+\tau s{{\displaystyle \sinh }}^{2}r}|}_{s=\tau =0}\mathrm{.}\end{array}

Comparing with the generating function of {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right), we obtain

\begin{array}{}\text{(42)}& N_m={\left(-\frac{1}{4}\sinh 2r\right)}^m{ℌ}_{m\mathrm{,}m}\left(\mathrm{0,0};-2\tanh r\right)\mathrm{.}\end{array}

Indeed, the normalized factor N_m is also reformed as the form of Legendre polynomial P_m(x) [24, 37] via computing high-order differential and the new formula of the polynomials P_m(x) [38]. Hence, by comparing with Eq. (42), we find that the relation between the polynomials {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right) and P_m(x) yields

\begin{array}{}\text{(43)}& {ℌ}_{m\mathrm{,}m}\left(\mathrm{0,0};-2\tanh r\right)=m!(\text{i}2\text{sech}r{)}^m{P}_m(\text{i}\sinh r)\mathrm{.}\end{array}

Specifically, m=0, {ℌ}_{\mathrm{0,0}}\left(\mathrm{0,0};-2\tanh r\right)=1, thus N_0=1 as expected, however m=1, {ℌ}_{\mathrm{1,1}}\left(\mathrm{0,0};-2\tanh r\right)=-2\tanh r, so N_1={{\displaystyle \sinh }}^{2}r, which is the normalized factor of the state aS\left(r\right)|0\rangle.

Another example is the photon-modulated state {\left(ta+r{a}^{†}\right)}^m|\gamma \rangle, which can be obtained via extending an elementary coherent superposition of photon subtraction and addition to the case of m order [39] and acting it on a coherent state |\gamma \rangle, and where the ratios t\mathrm{,}r are real and yield t^2+r^2=1. Indeed, this state can be realized physically via inputting the coherent state |\gamma \rangle into a parametric down amplifier and carrying out the detection of a single photon only at one of the output ports of the second beamsplitter in a optical experiment [39]. Using the normal ordering of {\left(ta+r{a}^{†}\right)}^m [24], we reform the state {\left(ta+r{a}^{†}\right)}^m|\gamma \rangle as

\begin{array}{}\text{(44)}& {\left(ta+r{a}^{†}\right)}^m|\gamma \rangle ={\left(-\text{i}\sqrt{\frac{rt}{2}}\right)}^m{H}_m\left(\text{i}\frac{t\gamma +r{a}^{†}}{\sqrt{2rt}}\right)|\gamma \rangle \mathrm{.}\end{array}

Thus, using the completeness of coherent states [36] and the generating function of {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right), we find that the normalized factor of {\left(ta+r{a}^{†}\right)}^m|\gamma \rangle reads

\begin{array}{}\text{(45)}& \begin{array}{ll}{\mathcal{N}}_m\hfill & ={\left(\frac{rt}{2}\right)}^m{\displaystyle \int }\frac{{\text{e}}^2\beta }{\pi }{\text{e}}^{-{\left|\beta \right|}^2-{\left|\gamma \right|}^2+{\beta }^{\ast }\gamma +\beta {\gamma }^{\ast }}\hfill \\ \hfill & \times H_m\left(-\text{i}\frac{t{\gamma }^{\ast }+r\beta }{\sqrt{2rt}}\right)H_m\left(\text{i}\frac{t\gamma +r{\beta }^{\ast }}{\sqrt{2rt}}\right)\hfill \\ \hfill & ={\left(\frac{rt}{2}\right)}^m{\frac{{\partial }^{2m}}{\partial s^m\partial {\tau }^m}{\text{e}}^{-s^2-{\tau }^2+2rs\tau /t+\varrho s+{\varrho }^{\ast }\tau }|}_{s=\tau =0}\hfill \\ \hfill & ={\left(\frac{rt}{2}\right)}^m{ℌ}_{m\mathrm{,}m}\left(\varrho \mathrm{,}{\varrho }^{\ast };\frac{2r}{t}\right)\mathrm{,}\hfill \end{array}\end{array}

where \varrho =\text{i}\sqrt{2}\left(t\gamma +r{\gamma }^{\ast }\right)/\sqrt{rt}. Particularly, when m=0, {ℌ}_{\mathrm{0,0}}\left(\varrho \mathrm{,}{\varrho }^{\ast };2r/t\right)=1, thus {\mathcal{N}}_0=1. For the case of m=1, {ℌ}_{\mathrm{1,1}}\left(\varrho \mathrm{,}{\varrho }^{\ast };2r/t\right)=2r/t+{\left|\varrho \right|}^2, so {\mathcal{N}}_1=r^2+{\left|t\gamma +r{\gamma }^{\ast }\right|}^2. On the other hand, when \gamma =0, {\mathcal{N}}_m becomes {\left(rt/2\right)}^m{ℌ}_{m\mathrm{,}m}\left(\mathrm{0,0};2r/t\right), which is exactly the normalized factor of the state {(ta+r{a}^{†})}^m|0\rangle. Indeed, this factor can be verified by using the newly derived integral formula (37) to give the operator identity as

\begin{array}{}\text{(46)}& \begin{array}{l}{H}_m\left(-\text{i}\frac{\sqrt{r}a}{\sqrt{2t}}\right)H_m\left(\text{i}\frac{\sqrt{r}{a}^{†}}{\sqrt{2t}}\right)\\ ={\displaystyle \int }\frac{{\text{d}}^2\alpha }{\pi }:{\text{e}}^{-\left(\alpha -a\right)\left({\alpha }^{\ast }-a^{†}\right)}\\ \times H_m\left(-\text{i}\frac{\sqrt{r}\alpha }{\sqrt{2t}}\right)H_m\left(\text{i}\frac{\sqrt{r}{\alpha }^{\ast }}{\sqrt{2t}}\right):\\ =:{ℌ}_{m\mathrm{,}m}\left(-\text{i}\frac{\sqrt{2r}}{\sqrt{t}}a\mathrm{,}\text{i}\frac{\sqrt{2r}}{\sqrt{t}}{a}^{†};\frac{2r}{t}\right):\mathrm{.}\end{array}\end{array}

Besides, for a two-mode quantum state H_{n\mathrm{,}m}(f{a}^{†}，g{b}^{†})S_2(r)|00\rangle, where H_{n\mathrm{,}m}\left(f{a}^{†}\mathrm{,}g{b}^{†}\right) are two-variable operator Hermite polynomials of orders(n, m), f, g are arbitrary real parameters and S_2(r)=\exp \left[r\left(a^{†}{b}^{†}-ab\right)\right] is two-mode squeezing operator [40, 41], we use the new special polynomials F_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon ) to address its normalization. Experimentally, this state can be generated by inputting the squeezed vacuum S_2(r)|00\rangle into two tunable paralleled beamsplitters and implementing the conditional multi-photon measurement in one of the output ports of each beamsplitter [42]. Here, its normalization factor reads

\begin{array}{}\text{(47)}& \begin{array}{ll}{\mathfrak{N}}_{n\mathrm{,}m}=\hfill & {\text{sech}}^{2}r\langle 00|{\text{e}}^{ab\tanh r}{H}_{n\mathrm{,}m}\left(fa\mathrm{,}gb\right)\hfill \\ \hfill & \times H_{n\mathrm{,}m}\left(f{a}^{†}\mathrm{,}g{b}^{†}\right){\text{e}}^{a^{†}{b}^{†}\tanh r}|00\rangle \mathrm{.}\hfill \end{array}\end{array}

Inserting the completeness of coherent states |\alpha \beta \rangle and using the integral formula (40), we get

\begin{array}{}\text{(48)}& \begin{array}{ll}{\mathfrak{N}}_{n\mathrm{,}m}=& \frac{{\partial }^{n+m}}{\partial s^n\partial {\tau }^m}\frac{{\partial }^{n+m}}{\partial {s^{\prime }}^n\partial {{\tau }^{\prime }}^m}\exp [-\frac{s\tau +s^{\prime }{\tau }^{\prime }}{\mathfrak{f}}\hfill \\ \hfill & {+f^{2}s{s}^{\prime }{{\displaystyle \cosh }}^{2}r+g^2\tau {\tau }^{\prime }{{\displaystyle \cosh }}^{2}r]|}_{s=\tau =s^{\prime }={\tau }^{\prime }=0}\mathrm{,}\hfill \end{array}\end{array}

where \mathfrak{f}=2/(2-fg\sinh 2r). Further, using the definition of the special polynomials {\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon ), thus the factor {\mathfrak{N}}_{n\mathrm{,}m} is simplified as

\begin{array}{}\text{(49)}& {\mathfrak{N}}_{n\mathrm{,}m}=\frac{{\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}n\mathrm{,}m}\left(\mathrm{0,0,0,0};\mathfrak{f}{f}^2{{\displaystyle \cosh }}^{2}r\mathrm{,}\mathfrak{f}{g}^2{{\displaystyle \cosh }}^{2}r\right)}{{\mathfrak{f}}^{n+m}}\mathrm{,}\end{array}

which is just proportional to a six-variable special polynomial. For r=0, \mathfrak{f}=\mathrm{1,}{\mathfrak{N}}_{n\mathrm{,}m} reduces to {\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}n\mathrm{,}m}\left(\mathrm{0,0,0,0};f^2\mathrm{,}{g}^2\right), corresponding to the normalized factor of H_{n\mathrm{,}m}\left(f{a}^{†}\mathrm{,}g{b}^{†}\right)|00\rangle, however for f=g=0, the factor {\mathfrak{N}}_{n\mathrm{,}m}={\left(-1\right)}^{n+m}n!m!.

4.2 Photocount distribution

The quantum photocount distribution was first introduced by Kelley and Kleiner [43], it refers to a probability statistical distribution of registering n photoelectrons within an interval. We here continue to use the multi-variable special polynomials to analytically investigate the photocount distribution of the state a^{m}S\left(r\right)|0\rangle and the evolution of the photocount distribution of the number state |m\rangle in the thermal channel.

For a single-mode quantum state ρ, the new formula for calculating its photocount distribution is defined as [44]

\begin{array}{}\text{(50)}& \begin{array}{ll}\mathcal{P}(n)=\hfill & \frac{{\xi }^n}{{\left(\xi -1\right)}^n}{\displaystyle \int }\frac{{\text{d}}^2\alpha }{\pi }{\text{e}}^{-\xi {\left|\alpha \right|}^2}{L}_n\left({\left|\alpha \right|}^2\right)\hfill \\ \hfill & \times \langle \sqrt{1-\xi }\alpha |\rho |\sqrt{1-\xi }\alpha \rangle \mathrm{,}\hfill \end{array}\end{array}

which is indeed related to the the average value of ρ under the coherent state |\sqrt{1-\xi }\alpha \rangle (phase-space Q distribution function), and where ξ is quantum detection probability of single photon in a certain interval. For the ideal case, i.e., \xi =1, \mathcal{P}(n) becomes the photon number distribution of the state ρ. Using the normal ordering of the density operator of the state a^{m}S\left(r\right)|0\rangle [37], we easily obtain the average value

\begin{array}{}\text{(51)}& \begin{array}{l}\langle \sqrt{1-\xi }\alpha |\rho |\sqrt{1-\xi }\alpha \rangle \\ =\frac{\text{sech}r{{\displaystyle \tanh }}^{m}r}{N_m{2}^m}{\text{e}}^{-\left(1-\xi \right){\left|\alpha \right|}^2+\frac{1-\xi }{2}\left({\alpha }^{\ast 2}+{\alpha }^2\right)\tanh r}\\ \times H_m\left(\varpi {\alpha }^{\ast }\right)H_m\left({\varpi }^{\ast }\alpha \right)\mathrm{,}\end{array}\end{array}

where \varpi =\text{i}{\left[\left(1-\xi \right)\tanh r\right]}^{\mathrm{1/2}}/2. Substituting Eq. (51) into Eq. (50) and using the relation {(-1)}^{n}n!L_n\left({\left|\alpha \right|}^2\right)=H_{n\mathrm{,}n}\left(\alpha \mathrm{,}{\alpha }^{\ast }\right), where L_n(\cdot ) is the Laguerre polynomial, and the generating functions of the one- and two-variable Hermite polynomials, we have

\begin{array}{}\text{(52)}& \begin{array}{ll}\mathcal{P}(n)=\hfill & \frac{{\xi }^n\text{sech}r{{\displaystyle \tanh }}^{m}r}{N_m{2}^{m}n!{\left(1-\xi \right)}^n}\frac{{\partial }^{2m}}{\partial s^m\partial {\tau }^m}\frac{{\partial }^{2n}}{\partial {s^{\prime }}^n\partial {{\tau }^{\prime }}^n}\hfill \\ \hfill & \times {\text{e}}^{-s^2-{\tau }^2-s^{\prime }{\tau }^{\prime }}{\displaystyle \int }\frac{{\text{d}}^2\alpha }{\pi }{\text{e}}^{-{\left|\alpha \right|}^2+\frac{1-\xi }{2}\left({\alpha }^{\ast 2}+{\alpha }^2\right)\tanh r}\hfill \\ \hfill & {\times {\text{e}}^{2s\varpi {\alpha }^{\ast }+2\tau {\varpi }^{\ast }\alpha +s^{\prime }\alpha +{\tau }^{\prime }{\alpha }^{\ast }}|}_{s=\tau =s^{\prime }={\tau }^{\prime }=0}\mathrm{.}\hfill \end{array}\end{array}

Further, using the integral formula (30) in Ref. [37], we obtain

\begin{array}{}\text{(53)}& \begin{array}{ll}\mathcal{P}(n)=\hfill & \frac{{\xi }^n{\omega }^{\mathrm{1/2}}\text{sech}r{{\displaystyle \tanh }}^{m}r}{N_m{2}^{m}n!{\left(1-\xi \right)}^n}\frac{{\partial }^{2m}}{\partial s^m\partial {\tau }^m}\frac{{\partial }^{2n}}{\partial {s^{\prime }}^n\partial {{\tau }^{\prime }}^n}\hfill \\ \hfill & \exp [\left(\omega -1\right)s^{\prime }{\tau }^{\prime }+4\omega {\left|\varpi \right|}^{2}s\tau +2\omega {\varpi }^{\ast }\tau {\tau }^{\prime }\hfill \\ \hfill & +2\omega \varpi s{s}^{\prime }-4\epsilon \omega \varpi s{\tau }^{\prime }-4\epsilon \omega {\varpi }^{\ast }{s}^{\prime }\tau -\epsilon \omega {s^{\prime }}^2\hfill \\ \hfill & -\epsilon \omega {{\tau }^{\prime }}^2-\left(1+4\epsilon \omega {\varpi }^2\right)s^2\hfill \\ \hfill & -\left(1+4\epsilon \omega {\varpi }^{\ast 2}\right){\tau }^2]|{}_{s=\tau =s^{\prime }={\tau }^{\prime }=0}\mathrm{,}\hfill \end{array}\end{array}

where we have set \omega =1/\left(1-4{\epsilon }^2\right) and \epsilon =\left[\left(\xi -1\right)\tanh r\right]/2. Noting the definition of {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right), we rewrite Eq. (53) as

\begin{array}{}\text{(54)}& \begin{array}{ll}\mathcal{P}(n)=\hfill & \frac{{\xi }^n{\omega }^{\mathrm{1/2}}\text{sech}r{{\displaystyle \tanh }}^{m}r}{N_m{2}^{m}n!{\left(1-\xi \right)}^n}\hfill \\ \hfill & \times {\displaystyle \sum _{l\mathrm{,}k\mathrm{,}{l}^{\prime }\mathrm{,}{k}^{\prime }=0}^{\infty }}\frac{{\varpi }^{k+l^{\prime }}{\varpi }^{\ast l+k^{\prime }}{\left(2\omega \right)}^{l+k}{\left(-4\omega \epsilon \right)}^{l^{\prime }+k^{\prime }}}{l!k!l^{\prime }!k^{\prime }!}\hfill \\ \hfill & \times \frac{{\partial }^{2m}}{\partial s^m\partial {\tau }^m}\frac{{\partial }^{2n}}{\partial {s^{\prime }}^n\partial {{\tau }^{\prime }}^n}{s}^{k+l^{\prime }}{\tau }^{l+k^{\prime }}{s}^{\text{'}k+k^{\prime }}{\tau }^{\text{'}l+l^{\prime }}\hfill \\ \hfill & \times \exp [\left(\omega -1\right)s^{\prime }{\tau }^{\prime }+4\omega {\left|\varpi \right|}^{2}s\tau -\epsilon \omega {s^{\prime }}^2-\epsilon \omega {{\tau }^{\prime }}^2\hfill \\ \hfill & -{\lambda }^{-2}{s}^2-{\lambda }^{\ast -2}{\tau }^2+sx+\tau y+s^{\prime }{x}^{\prime }\hfill \\ \hfill & {+y^{\prime }{\tau }^{\prime }]|}_{s=\tau =s^{\prime }={\tau }^{\prime }=x=y=x^{\prime }=y^{\prime }=0}\mathrm{,}\hfill \end{array}\end{array}

where \lambda =1/\sqrt{1+4{\varpi }^2\omega \epsilon } and we have added the exponential term {\text{exp}(sx+\tau y+s^{\prime }{x}^{\prime }+y^{\prime }{\tau }^{\prime })|}_{x=y=x^{\prime }=y^{\prime }=0} in order to meet the definition of {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right). Using the differential relation of {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right) in (14), we have

\begin{array}{}\text{(55)}& \begin{array}{ll}\mathcal{P}(n)=\hfill & \frac{{\xi }^n{\omega }^{\mathrm{1/2}}\text{sech}r{{\displaystyle \tanh }}^{m}r}{N_m{2}^m{\left(1-\xi \right)}^n{\left|\lambda \right|}^{2m}}n!{\left(m!\right)}^2{\left(\epsilon \omega \right)}^n\hfill \\ \hfill & \times {\displaystyle \sum _{l\mathrm{,}k\mathrm{,}{l}^{\prime }\mathrm{,}{k}^{\prime }=0}^{\infty }}\frac{{ϵ}^{k+l^{\prime }}{ϵ}^{\ast l+k^{\prime }}{\left(2\omega \right)}^{l+k}{\left(-4\omega \epsilon \right)}^{l^{\prime }+k^{\prime }}}{l!k!l^{\prime }!k^{\prime }!\left(m-l-k^{\prime }\right)!}\hfill \\ \hfill & \times \frac{{ℌ}_{m-k-l^{\prime }\mathrm{,}m-l-k^{\prime }}\left(\mathrm{0,0};4{\left|\lambda \varpi \right|}^2\omega \right)}{\left(m-k-l^{\prime }\right)!\left(n-l-l^{\prime }\right)!\left(n-k-k^{\prime }\right)!}\hfill \\ \hfill & \times {ℌ}_{n-k-k^{\prime }\mathrm{,}n-l-l^{\prime }}\left(\mathrm{0,0};\frac{\omega -1}{\epsilon \omega }\right)\mathrm{,}\hfill \end{array}\end{array}

which is the analytical photocount distribution of the state a^{m}S\left(r\right)|0\rangle, related to the product of two three-variable special polynomials, and where ϵ=\lambda \varpi /\sqrt{\epsilon \omega }. However, in the previous work the distribution \mathcal{P}(n) can only be given in the form of high-order partial differential that can’t be calculated [37]. When m=0, we obtain the photocount distribution of the state S\left(r\right)|0\rangle as

\begin{array}{}\text{(56)}& {\mathcal{P}}_{m=0}(n)=\frac{{\xi }^n{\omega }^{\mathrm{1/2}}{\left(\epsilon \omega \right)}^n\text{sech}r}{n!{\left(1-\xi \right)}^n}{ℌ}_{n\mathrm{,}n}\left(\mathrm{0,0};\frac{\omega -1}{\epsilon \omega }\right)\mathrm{,}\end{array}

which is consistent with the previous result after simple calculation [44]. For \xi =1, \epsilon =\varpi =0, \lambda =\omega =1, thus \mathcal{P}(n) reduces to the photon number distribution of the state a^{m}S\left(r\right)|0\rangle, as shown in Ref. [24].

Next, noting that the Wigner distribution evolution of the state |m\rangle for thermal noise reads

\begin{array}{}\text{(57)}& \begin{array}{ll}{W}_m\left(\alpha \mathrm{,}t\right)=\hfill & \frac{{\left[2\left(\overline{\mathfrak{n}}+1\right)\mathcal{T}-1\right]}^m}{\pi {\left(2\overline{\mathfrak{n}}\mathcal{T}+1\right)}^{m+1}}{L}_m\left({\left|\mathfrak{g}\alpha \right|}^2\right)\hfill \\ \hfill & \times \exp \left(-\frac{2{\left|\alpha \right|}^2}{2\overline{\mathfrak{n}}\mathcal{T}+1}\right)\mathrm{,}\hfill \end{array}\end{array}

where \mathfrak{g}=2{\text{e}}^{-\kappa t}/\sqrt{\left(2\overline{\mathfrak{n}}\mathcal{T}+1\right)\left[1-2\left(\overline{\mathfrak{n}}+1\right)\mathcal{T}\right]}, \mathcal{T}=1-{\text{e}}^{-2\kappa t}, \kappa is the decay rate and \overline{\mathfrak{n}} is the average thermal photon number of the environment, so we introduce another new formula for calculating the photocount distribution of the state ρ that is interrelated with the Wigner distribution W\left(\alpha \right), which is represented as [44]

\begin{array}{}\text{(58)}& \begin{array}{ll}\mathcal{P}(n)=\hfill & \frac{4{\left(-\xi \right)}^n}{{\left(2-\xi \right)}^{n+1}}{\displaystyle \int }{\text{d}}^2\alpha {\text{e}}^{-2\xi {\left|\alpha \right|}^2/(2-\xi )}\hfill \\ \hfill & \times L_n\left(\frac{4{\left|\alpha \right|}^2}{2-\xi }\right)W\left(\alpha \right)\mathrm{.}\hfill \end{array}\end{array}

Substituting Eq. (57) into Eq. (58) and using a similar way of obtaining Eq. (52), thus we give the analytical photocount distribution evolution of the state |m\rangle for thermal noise as

\begin{array}{}\text{(59)}& \begin{array}{ll}\mathcal{P}(n\mathrm{,}t)=\hfill & \frac{4{\xi }^n{\left[1-2\left(\overline{\mathfrak{n}}+1\right)\mathcal{T}\right]}^m}{n!m!{\left(2-\xi \right)}^{n+1}{\left(2\overline{\mathfrak{n}}\mathcal{T}+1\right)}^{m+1}}\hfill \\ \hfill & \times \frac{{\partial }^{2n}}{\partial s^n\partial {\tau }^n}\frac{{\partial }^{2m}}{\partial {s^{\prime }}^m\partial {{\tau }^{\prime }}^m}{\text{e}}^{-s\tau -s^{\prime }{\tau }^{\prime }}{\displaystyle \int }\frac{{\text{d}}^2\alpha }{\pi }\hfill \\ \hfill & {{\text{e}}^{-{\mathfrak{g}}^{\prime }{\left|\alpha \right|}^2+\frac{2\left(\alpha s+{\alpha }^{\ast }\tau \right)}{\sqrt{2-\xi }}+\mathfrak{g}\alpha {\tau }^{\prime }+\mathfrak{g}{\alpha }^{\ast }{s}^{\prime }}|}_{s=\tau =s^{\text{'}}={\tau }^{\text{'}}=0}\mathrm{,}\hfill \end{array}\end{array}

where we have set {\mathfrak{g}}^{\prime }=2\xi /(2-\xi )+2/(2\overline{\mathfrak{n}}\mathcal{T}+1). Further, using the special polynomials {\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,}y\mathrm{,}{x}^{\prime }\mathrm{,}{y}^{\prime };\nu \mathrm{,}\upsilon ), we rewrite the photocount distribution \mathcal{P}(n\mathrm{,}t) as

\begin{array}{}\text{(60)}& \begin{array}{ll}\mathcal{P}(n\mathrm{,}t)=\hfill & \frac{4{\xi }^n{\left[1-2\left(\overline{\mathfrak{n}}+1\right)\mathcal{T}\right]}^m{\left[{\mathfrak{g}}^{\prime }\left(2-\xi \right)-4\right]}^n}{n!m!{\left(2-\xi \right)}^{2n+1}{\left(2\overline{\mathfrak{n}}\mathcal{T}+1\right)}^{m+1}{{\mathfrak{g}}^{\prime }}^{n+m+1}}\hfill \\ \hfill & \times {\left({\mathfrak{g}}^{\prime }-{\mathfrak{g}}^2\right)}^m{\mathfrak{F}}_{n\mathrm{,}n\mathrm{,}m\mathrm{,}m}\left(\mathrm{0,0,0,0};\mathfrak{h}\mathrm{,}\mathfrak{h}\right)\mathrm{,}\hfill \end{array}\end{array}

where \mathfrak{h}=2\mathfrak{g}/\sqrt{\left({\mathfrak{g}}^{\prime }-{\mathfrak{g}}^2\right)\left[{\mathfrak{g}}^{\prime }\left(2-\xi \right)-4\right]}. In particular, for t\to \infty \mathrm{,}\mathcal{T}\to \mathrm{1,}\mathfrak{g}\to \mathrm{0,}{\mathfrak{g}}^{\prime }\to 2\xi /(2-\xi )+2/(2\overline{\mathfrak{n}}+1), and \mathfrak{h}\to 0, thus \mathcal{P}(n\mathrm{,}t) reduces to the photocount distribution of the thermal field with the average photon number \overline{\mathfrak{n}}, that is, \mathcal{P}(n\mathrm{,}\infty )\to {\left(\xi \overline{\mathfrak{n}}\right)}^n/{\left(\xi \overline{\mathfrak{n}}+1\right)}^{n+1}, which agrees with the existed result in Ref. [44]. when \xi =1, the photocount distribution \mathcal{P}(n\mathrm{,}t) becomes the evolution of photon number distribution of the state |m\rangle for thermal noise, that is

\begin{array}{}\text{(61)}& \begin{array}{ll}{\mathcal{P}}_{\xi =1}(n\mathrm{,}t)=\hfill & \frac{4{\left[1-2\left(\overline{\mathfrak{n}}+1\right)\mathcal{T}\right]}^m{\left({\mathfrak{g}}^{″}-4\right)}^n{\left({\mathfrak{g}}^{″}-{\mathfrak{g}}^2\right)}^m}{n!m!{\left(2\overline{\mathfrak{n}}\mathcal{T}+1\right)}^{m+1}{{\mathfrak{g}}^{″}}^{n+m+1}}\hfill \\ \hfill & \times {\mathfrak{F}}_{n\mathrm{,}n\mathrm{,}m\mathrm{,}m}\left(\mathrm{0,0,0,0};{\mathfrak{h}}^{\prime }\mathrm{,}{\mathfrak{h}}^{\prime }\right)\mathrm{,}\hfill \end{array}\end{array}

where {\mathfrak{g}}^{″}=2+2/(2\overline{\mathfrak{n}}\mathcal{T}+1) and {\mathfrak{h}}^{\prime }=2\mathfrak{g}/\left[\left({\mathfrak{g}}^{″}-{\mathfrak{g}}^2\right)\left({\mathfrak{g}}^{″}-4\right)\right]. Further, letting t\to \infty in (61), {\mathfrak{g}}^{″}\to 2+2/(2\overline{\mathfrak{n}}+1), {\mathfrak{h}}^{\text{'}}\to 0, thus {\mathcal{P}}_{\xi =1}(n\mathrm{,}\infty )\to {\overline{\mathfrak{n}}}^n/{\left(\overline{\mathfrak{n}}+1\right)}^{n+1}, a Bose–Einstein statistical distribution.

4.3 Wigner distribution

Wigner distribution is a very useful tool in comprehensively describing nonclassical states from the phase-space point of view. Using the multi-variable special polynomials, we can simplify the calculations of the Wigner distributions of quantum states and conveniently investigate their nonclassicality in the phase space. Generally, for the state \rho, its Wigner distribution is defined as W(\alpha )=\text{tr}[\rho \text{Δ}(\alpha )] [45], where \text{Δ}(\alpha ) is the Wigner operator in the coherent state representation. Hence, using the identity(44), the state {\left(ta+r{a}^{†}\right)}^m|\gamma \rangle possesses the Wigner distribution

\begin{array}{}\text{(62)}& \begin{array}{l}W(\alpha )\\ =\frac{{\left(rt\right)}^m{\text{e}}^{2{\left|\alpha \right|}^2-{\left|\gamma \right|}^2}}{2^m{\mathcal{N}}_m}{\displaystyle \int }\frac{{\text{d}}^2{\alpha }^{\prime }}{{\pi }^2}{H}_m\left(-\text{i}\frac{t{\gamma }^{\ast }+r{\alpha }^{\prime }}{\sqrt{2rt}}\right)\\ \times H_m\left(\text{i}\frac{t\gamma -r{{\alpha }^{\prime }}^{\ast }}{\sqrt{2rt}}\right){\text{e}}^{-{\left|{\alpha }^{\prime }\right|}^2+\left(2\alpha -\gamma \right){{\alpha }^{\prime }}^{\ast }-\left(2{\alpha }^{\ast }-{\gamma }^{\ast }\right){\alpha }^{\prime }}\mathrm{.}\end{array}\end{array}

Further, using the generating function of H_n(x) and the integral formula (30) in Ref. [37], we have

\begin{array}{}\text{(63)}& \begin{array}{ll}W(\alpha )\hfill & =\frac{{\left(rt\right)}^m}{\pi 2^m{\mathcal{N}}_m}{\text{e}}^{-2{\left|\alpha -\gamma \right|}^2}\frac{{\partial }^{2m}}{\partial s^m\partial {\tau }^m}\exp (-s^2\hfill \\ \hfill & {-{\tau }^2-\frac{2r}{t}s\tau +\varkappa s+{\varkappa }^{\ast }\tau )|}_{s=\tau =0}\hfill \\ \hfill & =\frac{{\left(rt\right)}^m{\text{e}}^{-2{\left|\alpha -\gamma \right|}^2}}{\pi 2^m{\mathcal{N}}_m}{ℌ}_{m\mathrm{,}m}\left(\varkappa \mathrm{,}{\varkappa }^{\ast };-\frac{2r}{t}\right)\mathrm{,}\hfill \end{array}\end{array}

where \varkappa =\text{i}\sqrt{2}\left(t\gamma -r{\gamma }^{\ast }+2r{\alpha }^{\ast }\right)/\sqrt{rt}. Especially, when m=0, {ℌ}_{\mathrm{0,0}}\left(\varkappa \mathrm{,}{\varkappa }^{\ast };-2r/t\right)=1, so W_0(\alpha )={\pi }^{-1}{\text{e}}^{-2{\left|\alpha -\gamma \right|}^2}, which is the Wigner distribution of the coherent state |\gamma \rangle. However, for m=\mathrm{1,}{ℌ}_{11}\left(\varkappa \mathrm{,}{\varkappa }^{\ast };-2r/t\right)={\left|\varkappa \right|}^2-2r/t, thus we have

\begin{array}{}\text{(64)}& W_1(\alpha )=\frac{{\left|t\gamma -r{\gamma }^{\ast }+2r{\alpha }^{\ast }\right|}^2-r^2}{r^2+{\left|t\gamma +r{\gamma }^{\ast }\right|}^2}{W}_0(\alpha )\mathrm{,}\end{array}

which is in good agreement with the result in Ref. [39].

To clearly reveal how the coherent superposed operation {\left(ta+r{a}^{†}\right)}^m produces the nonclassicality of the coherent state |\gamma \rangle, we define the special polynomials {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right) as a new kind of “Wolfram” functions used like Hermite polynomials in the Mathematica programs, the result of which makes the writing easier and expedites calculation speed, so we more quickly simulate the Wigner distribution W(\alpha ). In Fig. 1, we plot the Wigner distribution W(\alpha ) for different parameters m\mathrm{,}r and \gamma, which shows that the negative volume of W(\alpha ) as a nonclassicality measure increases monotonically with the ratio r of the photon addition a^{†} and always decreases with the amplitude \gamma of the coherent state |\gamma \rangle, but constantly fluctuates with the order m. Hence, in order to enhance the nonclassicality of the state {\left(ta+r{a}^{†}\right)}^m|\gamma \rangle, we require a coherent superposition with an exact order m and a higher ratio r, and an initial coherent state with a smaller amplitude γ.

For the state H_{n\mathrm{,}m}\left(f{a}^{†}\mathrm{,}g{b}^{†}\right)S_2(r)|00\rangle, substituting its density operator and the coherent state representation of the Wigner operator \text{Δ}(\alpha \mathrm{,}\beta ) into the definition of the Wigner distribution, thus we obtain

\begin{array}{}\text{(65)}& \begin{array}{ll}W(\alpha \mathrm{,}\beta )=\hfill & \frac{{\text{sech}}^{2}r}{{\pi }^2{\mathfrak{N}}_{n\mathrm{,}m}}{\text{e}}^{2\left({\left|\alpha \right|}^2+{\left|\beta \right|}^2\right)}{\displaystyle \int }\frac{{\text{d}}^2{\alpha }^{\prime }{\text{d}}^2{\beta }^{\prime }}{{\pi }^2}\hfill \\ \hfill & H_{n\mathrm{,}m}\left(-f{{\alpha }^{\prime }}^{\ast }\mathrm{,}-g{{\beta }^{\prime }}^{\ast }\right)H_{n\mathrm{,}m}\left(f{\alpha }^{\prime }\mathrm{,}g{\beta }^{\prime }\right)\hfill \\ \hfill & \times \exp [\left({\alpha }^{\prime }{\beta }^{\prime }+{{\alpha }^{\prime }}^{\ast }{{\beta }^{\prime }}^{\ast }\right)\tanh r\hfill \\ \hfill & +2\left(\alpha {{\alpha }^{\prime }}^{\ast }+\beta {{\beta }^{\prime }}^{\ast }-{\left|{\alpha }^{\prime }\right|}^2\right)\hfill \\ \hfill & -2\left({\left|{\beta }^{\prime }\right|}^2+{\beta }^{\ast }{\beta }^{\prime }+{\alpha }^{\ast }{\alpha }^{\prime }\right)]\mathrm{.}\hfill \end{array}\end{array}

Using the generating function of H_{n\mathrm{,}m}(x\mathrm{,}y) and the formula (40) to integrate over {\alpha }^{\prime } and {\beta }^{\prime } in Eq. (65), we get the analytical Wigner distribution W(\alpha \mathrm{,}\beta ) as

\begin{array}{}\text{(66)}& \begin{array}{l}W(\alpha \mathrm{,}\beta )\\ =\frac{1}{{\pi }^2{\mathfrak{N}}_{n\mathrm{,}m}}{\text{e}}^{2\left({\left|\beta \right|}^2-{\left|\alpha \right|}^2\right)-{\left|\mathfrak{G}\right|}^2{\cosh }^{2}r}\frac{{\partial }^{n+m}}{\partial s^n\partial {\tau }^m}\frac{{\partial }^{n+m}}{\partial {s^{\prime }}^n\partial {{\tau }^{\prime }}^m}\\ \exp [-{\mathfrak{f}}^{-1}\left(s\tau +s^{\prime }{\tau }^{\prime }\right)-g^2\tau {\tau }^{\prime }{\cosh }^{2}r\\ -f^{2}s{s}^{\prime }{\cosh }^{2}r-g{\mathfrak{G}}^{\ast }\tau {\cosh }^{2}r-g\mathfrak{G}{\tau }^{\prime }{\cosh }^{2}r\\ +2f\cosh r\left(\alpha \cosh r-{\beta }^{\ast }\sinh r\right)s\\ +2f{\cosh r\left({\alpha }^{\ast }\cosh r-\beta \sinh r\right)s^{\prime }]|}_{s=\tau =s^{\prime }={\tau }^{\prime }=0}\\ =\frac{1}{{\pi }^2{\mathfrak{N}}_{n\mathrm{,}m}{\mathfrak{f}}^{n+m}}\exp [-2\left({\left|\alpha \right|}^2+{\left|\beta \right|}^2\right)\cosh 2r\mathrm{.}\\ +2\left(\alpha \beta +{\alpha }^{\ast }{\beta }^{\ast }\right)\sinh 2r]\\ \times {\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}n\mathrm{,}m}\left({\varrho }^{\prime }\mathrm{,}{\varkappa }^{\prime }\mathrm{,}{{\varrho }^{\prime }}^{\ast }\mathrm{,}{{\varkappa }^{\prime }}^{\ast };-f^2\mathfrak{f}{\cosh }^{2}r\mathrm{,}-g^2\mathfrak{f}{\cosh }^{2}r\right)\mathrm{,}\end{array}\end{array}

where we have set {\varrho }^{\prime }=2f\sqrt{\mathfrak{f}}\cosh r\left(\alpha \cosh r-{\beta }^{\ast }\sinh r\right)\mathrm{,}{\varkappa }^{\prime }=-g\sqrt{\mathfrak{f}}{\mathfrak{G}}^{\ast }{{\displaystyle \cosh }}^{2}r and \mathfrak{G}=2\alpha \tanh r-2{\beta }^{\ast }. Especially, for r=\mathrm{0,}\mathfrak{f}=\mathrm{1,}\mathfrak{G}=-2{\beta }^{\ast }\mathrm{,}{\varrho }^{\prime }=2f\alpha, and {\varkappa }^{\prime }=2g\beta, thus W(\alpha \mathrm{,}\beta ) becomes the Wigner distribution of the state H_{n\mathrm{,}m}\left(f{a}^{†}\mathrm{,}g{b}^{†}\right)|00\rangle, i.e.,

\begin{array}{}\text{(67)}& \begin{array}{ll}{W}_0(\alpha \mathrm{,}\beta )=\hfill & \frac{{\text{e}}^{-2\left({\left|\alpha \right|}^2+{\left|\beta \right|}^2\right)}}{{\pi }^2{\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}n\mathrm{,}m}\left(\mathrm{0,0,0,0};f^2\mathrm{,}{g}^2\right)}\hfill \\ \hfill & \times {\mathfrak{F}}_{n\mathrm{,}m\mathrm{,}n\mathrm{,}m}(2f\alpha \mathrm{,2}g\beta \mathrm{,2}f{\alpha }^{\ast }\mathrm{,}\hfill \\ \hfill & 2g{\beta }^{\ast };-f^2\mathrm{,}-g^2)\mathrm{.}\hfill \end{array}\end{array}

When f=g=0, {\varrho }^{\prime }={\varkappa }^{\prime }=0 and \mathfrak{f}=1, so W(\alpha \mathrm{,}\beta ) reduces to

\begin{array}{}\text{(68)}& \begin{array}{ll}{W}_{\mathrm{0,0}}(\alpha \mathrm{,}\beta )=\hfill & \frac{1}{{\pi }^2}\exp [-2\left({\left|\alpha \right|}^2+{\left|\beta \right|}^2\right)\cosh 2r\hfill \\ \hfill & +2\left(\alpha \beta +{\alpha }^{\ast }{\beta }^{\ast }\right)\sinh 2r]\mathrm{,}\hfill \end{array}\end{array}

which is just the Wigner distribution of two-mode squeezed vacuum [23].

Similarly, defining the special polynomials F_{n\mathrm{,}m\mathrm{,}l\mathrm{,}k}(x\mathrm{,} y\mathrm{,}{x}^{\text{'}}\mathrm{,}{y}^{\text{'}};\nu \mathrm{,}\upsilon ) as another new kind of “Wolfram” functions used in the Mathematica programs, we easily proceed numerical simulation to the Wigner distribution W(\alpha \mathrm{,}\beta ) and plot how the different parameters f\mathrm{,}g\mathrm{,}n\mathrm{,}mand r affect the nonclassicality of the state H_{n\mathrm{,}m}\left(f{a}^{†}\mathrm{,}g{b}^{†}\right)S_2(r)|00\rangle in Fig. 2. It is found that the distribution W(\alpha \mathrm{,}\beta ) exhibits two negative dips on both sides of the upward main peak for m=n\ne 0 and the negative dips gradually increase with the increase of m=n. However, the negativity shows the uncertain change of fluctuating with increasing m\ne n, which is also valid for the parameters f\mathrm{,}g. Besides, the Wigner distribution W(\alpha \mathrm{,}\beta ) is squeezed in a certain direction as a signal of nonclassicality and the negative region increases with the squeezing r. In sum, the increase of the parameters m=n and r can monotonically enhance the nonclassicality, but the increase of f\mathrm{,}g and m\ne n has irregular effects on it.

5 Conclusions

In this article, by introducing some commutative superposition operators and using the operator ordering method, we have introduced two new multi-variable special polynomials and their generating functions, and presented several new operator identities and integral formulas involving the two multi-variable special polynomials. It is found that, instead of calculating the high-order partial differential, the normalization, photocount distributions and Wigner distributions of several quantum states that can be generated physically can be simplified as the forms of multi-variable special polynomials, which bring us convenience for further studying their properties and applications. Especially, the relation between three-variable special polynomial {ℌ}_{n\mathrm{,}m}\left(x\mathrm{,}y;\vartheta \right) and Legendre polynomial P_m(x) is found, and the analytical photocount distribution of the state a^{m}S\left(r\right)|0\rangle is presented rather than high-order partial differential that cannot be calculated. Hence, we believe that the multi-variable special polynomials and other relevant identities developed in this article could be widely used in physics and mathematics theory like the single- and two-mode Hermite polynomials. Besides, we hope to further generalize the polynomials to the cases with more variables and indexes in the future work.

References

[1] Hoffmann E., Hussin V., Marquette I., and Zhang Y. Z., Nonclassical behaviour of coherent states for systems constructed using exceptional orthogonal polynomials, J. Phys. A Math. Theor. 51(8), 085202 (2018)

[2] Hwang W. and Ryoo C. S., Differential equations associated with two variable degenerate Hermite polynomial, Mathematics 8(2), 228 (2020)

[3] Fan Y. and Fan Y., New eigenmodes of propagation in quadratic graded index media and complex fractional fourier transform, Commum. Theor. Phys. 39(1), 97 (2003)

[4] Fan Y. and Klauder J. R., Weyl correspondence and Prepresentation as operator Fredholm equations and their solutions, J. Phys. Math. Gen. 39(34), 10849 (2006)

[5] Fan Y. and Chen J., Atomic coherent states studied by virtue of the EPR entangled state and their Wigner functions, Eur. Phys. J. D 23(3), 437 (2003)

[6] Meng G., Wang J. S., and Fan H. Y., Atomic coherent state as the eigenstates of a two-dimensional anisotropic harmonic oscillator in a uniform magnetic field, Mod. Phys. Lett. A 24(38), 3129 (2009)

[7] Fan Y., Wan Z. L., Wu Z., and Zhang P. F., A new kind of special function and its application, Chin. Phys. B 24(10), 100302 (2015)

[8] Wang S., Meng X. G., and Fan H. Y., Time evolution of angular momentum coherent state derived by virtue of entangled state representation and a new binomial theorem, Chin. Phys. B 28(10), 100301 (2019)

[9] Meng G., Liu J. M., Wang J. S., and Fan H. Y., New generalized binomial theorems involving two-variable Hermite polynomials via quantum optics approach and their applications, Eur. Phys. J. D 73(2), 32 (2019)

[10] Fan Y. and Jiang T. F., Two-variable Hermite polynomials as time-evolutional transition amplitude for driven harmonic oscillator, Mod. Phys. Lett. B 21(08), 475 (2007)

[11] Fan Y. and Fan Y., New eigenmodes of propagation in quadratic graded index media and complex fractional Fourier transform, Commun. Theor. Phys. 39(1), 97 (2003)

[12] Fan Y. and Xu X. F., Talbot effect in a quadraticindex medium studied with two-variable Hermite polynomials and entangled states, Opt. Lett. 29(10), 1048 (2004)

[13] Dodonov V., Asymptotic formulae for two-variable Hermite polynomials, J. Phys. Math. Gen. 27(18), 6191 (1994)

[14] Fan Y. and Fan Y., New bosonic operator ordering identities gained by the entangled state representation and two-variable Hermite polynomials, Commum. Theor. Phys. 38(3), 297 (2002)

[15] Meng G., Wang J. S., Yang Z. S., Zhang X. Y., Zhang Z. T., Liang B. L., and Li K. C., Squeezed Hermite polynomial state: Nonclassical features and decoherence behavior, J. Opt. 22(1), 015201 (2020)

[16] Fan Y. and Ye X., Hermite polynomial states in two-mode Fock space, Phys. Lett. A 175(6), 387 (1993)

[17] Bergou A., Hillery M., and Yu D. Q., Minimum uncertainty states for amplitude-squared squeezing: Hermite polynomial states, Phys. Rev. A 43(1), 515 (1991)

[18] Zhang L., Yuan H. C., Hu L. Y., and Xu X. X., Synthesis of Hermite polynomial excited squeezed vacuum states from two separate single-mode squeezed vacuum states, Opt. Commun. 356(23), 223 (2015)

[19] Górska K., Horzela A., and Szafraniec F. H., Holomorphic Hermite polynomials in two variables, J. Math. Anal. Appl. 470(2), 750 (2019)

[20] Casper W. R., Kolb S., and Yakimov M., Bivariate continuous q-Hermite polynomials and deformed quantum Serre relations, arXiv: 2002.07895 (2020)

[21] Fan Y., Lu H. L., and Fan Y., Newton–Leibniz integration for ket–bra operators in quantum mechanics and derivation of entangled state representations, Ann. Phys. 321(2), 480 (2006)

[22] Fan Y., Entanged states, squeezed states gained via the route of developing Dirac’s sybmolic method and their applications, Int. J. Mod. Phys. B 18 (10n11), 1387 (2004)

[23] Meng G., Wang J. S., Liang B. L., and Han C. X., Evolution of a two-mode squeezed vacuum for amplitude decay via continuous-variable entangled state approach, Front. Phys. 13(5), 130322 (2018)

[24] Meng G., Wang Z., Fan H. Y., and Wang J. S., Nonclassicality and decoherence of photon-subtracted squeezed vacuum states, J. Opt. Soc. Am. B 29(11), 3141 (2012)

[25] Meng G., Wang J. S., Liang B. L., and Du C. X., Optical tomograms of multiple-photon-added Gaussian states via the intermediate state representation theory, J. Exp. Theor. Phys. 127(3), 383 (2018)

[26] Fan Y. and Klauder J. R., Eigenvectors of two particles’ relative position and total momentum, Phys. Rev. A 49(2), 704 (1994)

[27] Fan Y., Zaidi H. R., and Klauder J. R., New approach for calculating the normally ordered form of squeeze operators, Phys. Rev. D 35(6), 1831 (1987)

[28] Wang S., Meng X. G., and Fan H. Y., s-parameterized Weyl transformation and the corresponding quantization scheme, Chin. Phys. B 24(1), 014203 (2015)

[29] Du X., Meng X. G., Zhang R., and Wang J. S., Analytical and numerical investigations of displaced thermal state evolutions in a laser process, Chin. Phys. B 26(12), 120301 (2017)

[30] Liu H., Zhang Y. B., Yu Y. F., and Zhang Z. M., Photon–phonon squeezing and entanglement in a cavity optomechanical system with a flying atom, Front. Phys. 14(1), 12601 (2019)

[31] Liu J., Ye W., Zhou W. D., Zhang H. L., Huang J. H., and Hu L. Y., Entanglement of coherent superposition of photon-subtraction squeezed vacuum, Front. Phys. 12(5), 120307 (2017)

[32] Dodonov V. V., “Nonclassical” states in quantum optics: A “squeezed” review of the first 75 years, J. Opt. B: Quantum Semiclass. Opt. 4(1), R1 (2002)

[33] Wenger J., Tualle-Brouri R., and Grangier P., Non-Gaussian statistics from individual pulses of squeezed light, Phys. Rev. Lett. 92(15), 153601 (2004)

[34] Wakui K., Takahashi H., Furusawa A., and Sasaki M., Photon subtracted squeezed states generated with periodically poled KTiOPO₄, Opt. Express 15(6), 3568 (2007)

[35] Fan Y., New antinormal ordering expansion for density operators, Phys. Lett. A 161(1), 1 (1991)

[36] Glauber R., Coherent and incoherent states of the radiation field, Phys. Rev. 131(6), 2766 (1963)

[37] Hu Y. and Fan H. Y., Statistical properties of photonsubtracted squeezed vacuum in thermal environment, J. Opt. Soc. Am. B 25(12), 1955 (2008)

[38] Fan Y., Meng X. G., and Wang J. S., New form of Legendre polynomials obtained by virtue of excited squeezed state and IWOP technique in quantum optics, Commum. Theor. Phys. 46(5), 845 (2006)

[39] Lee Y. and Nha H., Quantum state engineering by a coherent superposition of photon subtraction and addition, Phys. Rev. 82(5), 053812 (2010)

[40] Zhong R., Wang X., and Qin W., Towards quantum entanglement of micromirrors via a two-level atom and radiation pressure, Front. Phys. 13(5), 130319 (2018)

[41] Liu T., Zheng Z. F., Zhang Y., Fang Y. L., and Yang C. P., Transferring entangled states of photonic cat-state qubits in circuit QED, Front. Phys. 15(2), 21603 (2020)

[42] Yuan C., Xu X. X., and Xu Y. J., Generating two variable Hermite polynomial excited squeezed vacuum states by conditional measurement on beam splitters, Optik 172(21), 1034 (2018)

[43] Kelley L. and Kleiner W. H., Theory of electromagnetic field measurement and photoelectron counting, Phys. Rev. 136(2A), A316 (1964)

[44] Fan Y. and Hu L. Y., Two quantum-mechanical photocount formulas, Opt. Lett. 33(5), 443 (2008)

[45] Scully M. O. and Zubairy M. S., Quantum Optics, Cambridge: Cambridge University Press, 1997

Rights & permissions

Higher Education Press 2020