Multimode Bogoliubov transformation and Husimi's Q-function

In this paper, we present numerical schemes for evaluating the matrix elements of Gaussian/non-Gaussian operators in the Fock state basis, which are identified as multivariate Hermite polynomials (MHPs). Using the integral transformation operator to perform the multimode Bogoliubov transformation, Husimi's Q-functions of Gaussian/non-Gaussian operators are easily derived as the generating functions of MHPs.


Introduction
Multivariate Hermite polynomials (MHPs) play important roles in various field of research, including quantum optics [1,2] and molecular spectroscopy [3,4]. However, the evaluation of MHPs is notoriously difficult, and various algorithms have been developed for this task, such as Willink's recursion formula [5]. Willink [5] also showed that the MHPs could be converted to multivariate Gaussian moments (MGMs). Kan [6] developed an efficient direct summation formula for MGMs. In his paper, the connection between the MGMs and the hafnian [7,8,9,10] is also given. Therefore, the MHPs, MGMs, and the hafnian are mutually convertible quantities that have the same computational complexity, but are labeled differently.
In molecular vibronic spectroscopy, the vibronic transition amplitudes between the two different electronic potential energy surfaces (represented as harmonic potential wells) are called the Franck-Condon integral (FCI), which is usually evaluated as MHPs. Although all MHP parameters are real values, the computational difficulties have led to various numerical strategies [4,11] developed in computational chemistry packages. Doktorov et al. [12] interpreted the molecular vibronic transition as a Gaussian quantum optical transformation-the multimode Bogoliubov transformation. The authors decomposed the vibronic transition operation into elementary quantum optical operators (displacement, rotation, and squeezing). Proceeding from the work of Doktorov et al. [12], after 40 years, Huh et al. [13,14] connected the story of molecular vibronic spectroscopy to Gaussian boson sampling [15,7,8], which is a quantum sampling problem of Gaussian light. Some experiments were performed successfully by quantum devices at the small scale [16,17,18]. The properties of Gaussian states are well explained in Refs. [19,20], and the detailed calculations involving Gaussian/non-Gaussian operations can be found in Refs. [21,22,23]. Therein, Husimi's Q-function [24] for the Gaussian light was derived, as a generating function, from Wigner's function by convolution. Only a few papers [1] have shown the complex Bogoliubov parameter dependence on Husimi's Q-function explicitly.
In this paper, we connect the matrix elements of non-Gaussian operators in the Fock state basis to the MHPs so that the evaluation can be conducted in an equal footing as the Gaussian case. Unlike the recent papers [21,22,23], we start from the multimode Bogoliubov transformation to define the pure Gaussian operator. We adapt the integral operator method of Fan et al. [25,26] and follow the derivation steps in Ref. [25] throughout the paper to obtain Husimi's Q-function, which is used as a generating function of the MHPs for the Gaussian/non-Gaussian matrix elements in the Fock state basis. With the aid of the integral operator method proposed by Fan et al., the explicit parameter dependence on the multimode Bogoliubov transformation matrices is made in Husimi's Q-function.
Our work is presented as follows: The relationship between the matrix elements of Gaussian/non-Gaussian operators in the Fock state basis, the MHPs, and the multivariate Gaussian distribution is given first (some contents are taken from the PhD thesis (Ch. 3) of the author [11] for the following section). Then, the integral operator transformation [25,26] method is introduced and applied to the multimode Bogoliubov transformation for the corresponding Husimi's Q-function of Gaussian/non-Gaussian operators. Finally, we summarize the paper in the conclusion section.
2 Matrix elements of Gaussian/non-Gaussian operator, multivariate Hermite polynomials, and moments of multivariate Gaussian distribution The matrix elements of a multimode Gaussian operator (Ô G ) in the N -dimensional Fock states (|n = |n 1 , · · · , n N and |m = |m 1 , · · · , m N ) can be generated by Husimi's Q-function of a Gaussian operator (π −N α|Ô G |α ) [24]; this can easily be seen via the expansion of the coherent states (|α ) in the Fock state basis, as follows, with the unnormalized Husimi Q-function, Therefore, the matrix elements in the Fock basis are given as partial derivatives with respect to the phase parameters of the coherent states, i.e.
Here, the small and large letters in bold fonts are used to indicate column vectors and square matrices, respectively. The unnormalized Husimi Q-function of the Gaussian operator (exp(α † α) α|Ô G |α ) has a closed form of the multivariate Gaussian distribution function, whereᾱ T = (α T α † ). The matrix elements are given as the MHPs H, where the M -dimensional Hermite polynomials are defined as follows, whereṽ = M j=1 v j , with a complex symmetric matrix Λ having a symmetric positive definite real part. We can evaluate the MHPs using the following recursion relation, (6) which was derived by Willink [5]. Willink also found the recursion relation for the MHPs via the relation between the multivariate normal (Gaussian) moments E and the MHPs H, i.e.
where y is a random variable vector in a multivariate normal distribution, N (y m , Λ −1 ) with its mean vector y m and covariance matrix Λ −1 . Exploiting the Magnus series expansion for products of variables: where h = v/2 − l, Kan [6] developed an efficient algorithm to evaluate the MGM, here [ṽ/2] is the greatest integer less than or equal toṽ/2. An (alternative) iterative evaluation scheme for the MHPs [11], may be evaluated efficiently by identifying y m = −iµ and Λ = W −1 in Eqs. (8) and (9), i.e.
Then we can evaluate the matrix elements in Eq. (4) with the relation (10) involving the summation of the univariate Hermite polynomials: the summation over the s part in Eq. (10) is the iterative expression of the univariate Hermite polynomials. The matrix elements of the non-Gaussian operator, for example, the position operator (Q) or momentum operator (P), can also be evaluated with Eqs. (6) and (10) using the following Husimi Q-functions with an auxiliary generating function parameter vector (λ), because these are again complex Gaussian functions with the following identities, In the subsequent section, Husimi's Q-functions are presented as complex Gaussian functions via the integral operator method [25,26].

Multimode Bogoliubov integral transformation operator and Husimi's Q-function
The most general pure Gaussian operatorÔ G can be obtained as the product of the displacement and multimode squeezing operators [27,28], which can be further decomposed as two rotations and a single squeezing operators. The resulting multimode Bogoliubov transformation for the bosonic annihilation (â,b) and creation (â † ,b † ) operators is given as follows, The transformation matrix K satisfies the symplectic condition and it results in the following symplectic identities, where I N and O N are an N -dimensional identity and square zero matrices, respectively. We use, in this paper, the convention of Ma and Rhodes [27] for transformation of boson operator vectorx asÂxB ≡ (Âx 1B , . . . ,Âx NB ) T . O G can be decomposed into the elementary quantum optical operators as [28,14], via the singular value decomposition of transformation matrices S and R (Bloch-Messiah reduction [28]): where Σ is a real diagonal matrix of squeezing parameters, and U L and U R are N × N unitary matrices. The elementary quantum optical unitary operators are defined as follows [27]: squeezing operator: These optical operators act onâ as, respectively: Fan et al. [25,26] proposed an integral form of the multimode squeezing operator, where d 2 z = N k dRe(z k )dIm(z k ), such thatÔ † SâÔ S = Sâ − Râ † .Ô G can also be given in an integral form by applying the displacement operatorD(t) to the integral operatorÔ S , i.e. where As shown by Fan et al. [25,26], the normal ordering ofÔ G is derived by using the integral form in Eq. (27) to find Husimi's Q-function. The coherent states in the integral operator are expressed as a vacuum projection form to exploit the normal ordering of the vacuum projection operator (|0 0|) for the integration, where ξ T = (z T z * T ), and with the symmetric matrix W, By introducing the vacuum projection operator identity [29,26], where the symbol :X: denotes the normal ordering of the operatorX, Eq. (27) becomes integrable, The Gaussian integral in Eq. (33), with a complex symmetric covariance matrix, can be evaluated using the following formula [26,30], where A and D are N × N symmetric matrices. After integration, we find a normal ordering of the Gaussian operator, which is again a Gaussian function of the bosonic operators: where we have used the symplectic identities (18) for Eq. (35) is further rearranged aŝ Finally, the unnormalized Husimi Q-function of the Gaussian operator is given as a complex Gaussian function, Eq. (40) agrees with the expression of Dondonov and Man'ko [1], which was obtained by a complex linear canonical transformation [31,32]. Alternatively, we can obtain the unnormalrized Husimi Q-function via the integral transformation directly without invoking the normally ordered form, that is, exp(α † α) α|Ô G |α = exp(α † α) |det(S)| d 2 z π N g(S, −R, t) α|Sz − Rz * + t z|α , We can exploit the direct integral transformation method to derive non-Gaussian matrix elements easily in the MHP formula. Here, we derive the matrix elements of the powers of position (Q = 1 √ 2 (â +â † )) and momentum (P = −i √ 2 (â −â † )) operators. Using the following identities [26], we express the exponential operators of position and momentum in the ordered form,