Bloch-Messiah reduction for twin beams of light

TL;DR

This paper analyzes the Bloch-Messiah reduction of pulsed parametric downconversion producing twin beams, revealing eigenvalue multiplicities, addressing mode ambiguity, and proposing methods for unique mode determination with practical examples.

Contribution

It introduces two approaches for uniquely determining squeezing eigenmodes in nondegenerate twin beams, addressing eigenvalue multiplicity and mode ambiguity issues.

Findings

Eigenvalues have at least double multiplicity in this regime

Modal functions can be derived from Schmidt modes or Hermitian eigenvalue problems

Eigenmode structure varies near phase-matching degeneracy

Abstract

We study the Bloch-Messiah reduction of parametric downconversion of light in the pulsed regime with a nondegenerate phase matching providing generation of twin beams. We find that in this case every squeezing eigenvalue has multiplicity at least two. We discuss the problem of ambiguity in the definition of the squeezing eigenmodes in this case and develop two approaches to unique determination of the latter. First, we show that the modal functions of the squeezing eigenmodes can be tailored from the Schmidt modes of the signal and idler beams. Alternatively, they can be found as a solution of an eigenvalue problem for an associated Hermitian squeezing matrix. We illustrate the developed theory by an example of frequency non-degenerate collinear twin beams generated in beta barium borate crystal. On this example we demonstrate how the squeezing eigenmodes can be approximated analytically on the basis of the Mehler's formula, extended to complex kernels. We show how the multiplicity of the eigenvalues and the structure of the eigenmodes are changed when the phase matching approaches the degeneracy in frequency.

Figures

Equations

$$A(\mathbf{r},t)=\sum _k{a}_k{f}_k^{\ast }(\mathbf{r},t)+\sum _k{a}_k^{†}{f}_k(\mathbf{r},t),$$

$$\mathbf{K}=\left(\begin{array}[]{cc}\mathbb{I}&0\\ 0&-\mathbb{I}\end{array}\right)$$

$${\mathbf{a}}^{\prime }=\mathbf{Sa},$$

$$\mathbf{S}=\left(\begin{array}[]{cc}\mathbb{S}_{0}&\mathbb{S}_{I}\\ \mathbb{S}_{I}^{*}&\mathbb{S}_{0}^{*}\end{array}\right),$$

$${\mathbf{a}}^{\prime }={\mathcal{U}}^{†}\mathbf{a}\mathcal{U}.$$

$$\mathbf{H}=\left(\begin{array}[]{cc}\mathbb{H}_{0}&\mathbb{H}_{I}\\ \mathbb{H}_{I}^{*}&\mathbb{H}_{0}^{*}\end{array}\right)$$

$$\mathbf{S}=\exp\left(\begin{array}[]{cc}-i\mathbb{H}_{0}&-i\mathbb{H}_{I}\\ i\mathbb{H}_{I}^{*}&i\mathbb{H}_{0}^{*}\end{array}\right).$$

$$\mathbf{S}=\left(\begin{array}[]{cc}e^{-i\mathbb{H}_{0}}&0\\ 0&e^{i\mathbb{H}_{0}^{*}}\end{array}\right)$$

$$\mathcal{H}$$

$$\mathcal{U}=e^{-i{\bar{a}}^{†}{\mathbb{H}}_0\bar{a}}.$$

$$\mathbf{S}=\exp\left(\begin{array}[]{cc}0&\mathbb{R}\\ \mathbb{R}&0\end{array}\right)=\left(\begin{array}[]{cc}\cosh(\mathbb{R})&\sinh(\mathbb{R})\\ \sinh(\mathbb{R})&\cosh(\mathbb{R})\end{array}\right).$$

$$\mathcal{H}$$

$$\mathcal{U}=e^{\frac{1}{2}\left({\bar{a}}^{†}\mathbb{R}{\bar{a}}^{†T}-{\bar{a}}^T\mathbb{R}\bar{a}\right)}.$$

$$\mathbf{S}=\left(\begin{array}[]{cc}\mathbb{V}&0\\ 0&\mathbb{V}^{*}\end{array}\right)\exp\left(\begin{array}[]{cc}0&\mathbb{R}\\ \mathbb{R}&0\end{array}\right)\left(\begin{array}[]{cc}\mathbb{Q}&0\\ 0&\mathbb{Q}^{*}\end{array}\right)^{\dagger},$$

$${\mathbb{S}}_0=\mathbb{V}\cosh (\mathbb{R}){\mathbb{Q}}^{†},{\mathbb{S}}_I=\mathbb{V}\sinh (\mathbb{R}){\mathbb{Q}}^T,$$

$$\mathcal{U}=e^{-i{\bar{a}}^{†}{\mathbb{H}}_V\bar{a}}{e}^{\frac{1}{2}\left({\bar{a}}^{†}\mathbb{R}{\bar{a}}^{†T}-{\bar{a}}^T\mathbb{R}\bar{a}\right)}{e}^{i{\bar{a}}^{†}{\mathbb{H}}_Q\bar{a}},$$

$$W(\mathbf{\alpha })=\frac{1}{{\pi }^{2n}}\int d^2\bar{\lambda }\mathrm{Tr}\left\{\rho e^{({\bar{a}}^{†}-{\bar{\alpha }}^{†})\bar{\lambda }-{\bar{\lambda }}^{†}(\bar{a}-\bar{\alpha })}\right\},$$

$$W(\mathbf{\alpha })=\frac{1}{(2\pi {)}^n\sqrt{\det \mathbf{\Sigma }}}{e}^{-\frac{1}{2}({\mathbf{\alpha }}^{†}-{\mathbf{\alpha }}_0^{†}){\mathbf{\Sigma }}^{-1}(\mathbf{\alpha }-{\mathbf{\alpha }}_0)},$$

$$\mathbf{\Sigma }=\frac{⟨(\mathbf{a}-{\mathbf{\alpha }}_0)({\mathbf{a}}^{†}-{\mathbf{\alpha }}_0^{†})⟩}{2}+\frac{⟨{\left[({\mathbf{a}}^{†}-{\mathbf{\alpha }}_0^{†}{)}^T(\mathbf{a}-{\mathbf{\alpha }}_0{)}^T\right]}^T⟩}{2}$$

$$\boldsymbol{\Sigma}=\left(\begin{array}[]{cc}\Sigma_{0}&\Sigma_{I}\\ \Sigma_{I}^{*}&\Sigma_{0}^{*}\end{array}\right),$$

$${\Sigma }_0=\frac{1}{2}⟨(\bar{a}-{\bar{\alpha }}_0)({\bar{a}}^{†}-{\bar{\alpha }}_0^{†})+{\left[({\bar{a}}^{†}-{\bar{\alpha }}_0^{†}{)}^T(\bar{a}-{\bar{\alpha }}_0{)}^T\right]}^T⟩$$

$${\Sigma }_I=⟨(\bar{a}-{\bar{\alpha }}_0)(\bar{a}-{\bar{\alpha }}_0{)}^T⟩$$

$$\mathbf{V}^{\dagger}\mathbf{a}^{\prime}=\left(\begin{array}[]{cc}\cosh(\mathbb{R})&\sinh(\mathbb{R})\\ \sinh(\mathbb{R})&\cosh(\mathbb{R})\end{array}\right)\mathbf{a}_{vac},$$

$$\mathbf{V}=\left(\begin{array}[]{cc}\mathbb{V}&0\\ 0&\mathbb{V}^{*}\end{array}\right),\,\,\,\mathbf{Q}=\left(\begin{array}[]{cc}\mathbb{Q}&0\\ 0&\mathbb{Q}^{*}\end{array}\right).$$

$$\mathcal{U}=e^{\frac{1}{2}\left({\bar{a}}^{†}\Gamma {\bar{a}}^{†T}-{\bar{a}}^T{\Gamma }^{\ast }\bar{a}\right)}{e}^{-i{\bar{a}}^{†}{\mathbb{H}}_V\bar{a}}{e}^{i{\bar{a}}^{†}{\mathbb{H}}_Q\bar{a}},$$

$$e^{-i{\bar{a}}^{†}{\mathbb{H}}_V\bar{a}}\bar{a}{e}^{i{\bar{a}}^{†}{\mathbb{H}}_V\bar{a}}=e^{i{\mathbb{H}}_V}\bar{a}={\mathbb{V}}^{†}\bar{a},$$

$${\mathbf{\Sigma }}^{\prime }=\frac{1}{2}\mathbf{S}{\mathbf{S}}^{†}=\frac{1}{2}\mathbf{V}{e}^{2\mathbf{R}}{\mathbf{V}}^{†},$$

$$\mathbf{R}=\left(\begin{array}[]{cc}0&\mathbb{R}\\ \mathbb{R}&0\end{array}\right).$$

$$E^{(+)}(t,z)=\frac{1}{2\pi }\int a(\Omega ,z)e^{-i\left({\omega }_0+\Omega \right)t}d\Omega ,$$

$$a(\Omega ,z)=ϵ(\Omega ,z)e^{ik(\Omega )\left(z-z_0\right)},$$

Full text

I Introduction

In the process of parametric downconversion (PDC) one photon of pump wave is converted into a pair of signal and idler photons, which are almost simultaneous Zel’dovich and Klyshko (1969); Burnham and Weinberg (1970). A multimode analysis Klyshko (1988); Hong and Mandel (1985) of this process shows that the correlation time between the two photons is determined by the inverse of the downconverted light bandwidth, which is typically in the sub-picosecond range. Thus, detection of the idler photon results in a localization of the signal one Zel’dovich and Klyshko (1969), which is known as “photon heralding” technique Hong and Mandel (1986), widely used today. Analysis of the joint quantum state of two photons shows that a discrete set of orthogonal modes, Schmidt modes, can be defined for each photon of the pair such that detection of a photon in an idler Schmidt mode projects the other photon onto the corresponding signal Schmidt mode Law et al. (2000); Grice et al. (2001). In the high-gain regime, when the initial spontaneous photon pairs experience amplification by stimulated emission and many photons are generated in the principal modes, nondegenerate PDC results in generation of two beams of light, signal and idler, having equal numbers of photons in any time interval longer than the inverse bandwidth of one beam. These highly correlated beams are known as “twin beams” and were generated both in the cavity Heidmann et al. (1987) and the single-pass Aytür and Kumar (1990) configurations. Early multimode theory of PDC Kolobov and Sokolov (1989); Reid (1989) allowed one to calculate many important properties of the generated field for the case of monochromatic pump. Later a modal decomposition for pulsed high-gain PDC was introduced Bennink and Boyd (2002), which represented the output field as an assembly of independent squeezing eigenmodes, each being in a single-mode squeezed state. Essentially similar approach was formulated in the language of symplectic transformations Braunstein (2005); McKinstrie and Karlsson (2013); Cariolaro and Pierobon (2016a, b), and is known as Bloch-Messiah reduction in reference to a similar reduction in the elementary particle physics Bloch and Messiah (1962); Balian et al. (1965). In the case of degenerate PDC, where the signal and the idler photons are indistinguishable, the Bloch-Messiah reduction has proven to be a powerful tool for determining the squeezing eigenmodes for single-pass optical parametric amplifiers (OPA) Wasilewski et al. (2006); Lvovsky et al. (2007) and multi-pass optical parametric oscillators Patera et al. (2010, 2012); Patera (2008); Jiang et al. (2012). Application of this formalism to a degenerate PDC with a monochromatic pump, undertaken recently by some of us Lipfert et al. (2018), resulted in a successful identification of bichromatic squeezing eigenmodes. For nondegenerate high-gain PDC it has been found Migdał and Wasilewski (2010) that up to certain level of gain the Schmidt modes of photon pairs determine a modal decomposition for the signal and idler beams such that the corresponding modes are in a two-mode squeezed state. This decomposition is obviously related to the Bloch-Messiah reduction but does not coincide with it, since in the latter, as introduced by Braunstein Braunstein (2005), the field is represented as an assembly of single-mode squeezed states.

The aim of the present article is to apply the Bloch-Messiah reduction to twin beams generated in a non-degenerate PDC with any pump, pulsed or continuous wave. We develop a general formalism applicable in all cases where signal and idler beams can be discriminated either by frequency, direction, or polarization. We find a fundamental result applicable to all these cases: every squeezing eigenvalue of twin beams has multiplicity at least two. This result follows directly from the symmetry of the interaction Hamiltonian in the case where the downconverted light is partitioned into a signal and an idler parts. As consequence, the squeezing eigenmodes are determined up to an orthogonal rotation in the subspace, corresponding to a given eigenvalue. Such a rotation leads to ambiguity in the definition of squeezing eigenmodes and often makes impossible obtaining reproducible results when the eigenmodes are found numerically by a singular value decomposition (SVD) of the Bogoliubov transformation matrices Braunstein (2005) or by a Takagi factorization of the transformation generator matrix Bennink and Boyd (2002); Cariolaro and Pierobon (2016a, b). We propose two solutions for this problem. First, we show how the squeezing eigenmodes can be composed from the Schmidt modes of each beam. Second, we show that the Takagi factorization of the squeezing matrix can be easily found from the spectral decomposition of an associated Hermitian matrix. As consequence, we find that under rather simple experimental conditions the modal functions of the squeezing eigenmodes can be identified as transform-limited waveforms. These modal functions create a basis in the subspace, corresponding to a given eigenvalue, and are represented (up to constant phase) by real functions of their spatio-temporal arguments. Measurement of pulsed squeezing is realized by pulsed homodying Medeiros de Araújo et al. (2014), which requires a precise shaping of the local oscillator pulse to fit the target eigenmode. Thus, finding a simple transform-limited basis in a continuous family of eigenfunctions means a significant simplification of the measurement process, especially in the full three-dimensional case, where the modal functions depend on two transversal coordinates and time. We illustrate our general results by an example where the modal functions depend on frequency only, postponing the treatment of the much more complicated three-dimensional case to a future publication.

The article is structured as follows. In Sec. II we give a brief review of the complex symplectic formalism for description of Gaussian unitary transformations, encompassing all variants of PDC with undepleted pump. The Bloch-Messiah reduction has a simple mathematical form in terms of this formalism. This formalism simplifies the mathematics of subsequent sections. Besides, we plan to apply this formalism in further publications to the description of three-dimensional spatio-temporal modes of a noncollinear OPA, which explains its rather detailed character. In Sec. III we review different ways of determining the squeezing eigenmodes and show their equivalence. We consider here the Magnus expansion of the field evolution operator. The complex symplectic formalism helps us to clarify the general structure of the lowest orders of the Magnus expansion and the role of the interaction picture reference frame. Section IV is devoted to generation of twin beams, which is a special case of PDC where the signal and the idler modes are distinct. We show that in this special case the Hamiltonian of the unitary transformation possesses a symmetry, leading to double multiplicity of the squeezing eigenvalues. We describe two procedures allowing one to avoid the ambiguity in the determination of the squeezing eigenmodes. The developed theory is illustrated in Sec. V by a concrete example of bright twin beams discriminated in frequency only with the conditions close to the experiment, which is underway in the Sorbonne University group. Section VI summarizes the results and makes an outlook for the future research.

II Symplectic formalism for Gaussian unitary transformations

II.1 Field decomposition

In quantum theory a multimode optical field is described by a Hermitian operator $A(\mathbf{r},t)$ of its vector potential in point $\mathbf{r}$ and time $t$. In the absence of sources this operator, in the Coulomb gauge, can be uniquely decomposed into a complete orthonormal set of modal functions $\{f_{k}(\mathbf{r},t)\}$ Grynberg et al. (2010); Fabre :

[TABLE]

so that $-\partial_{t}A(\mathbf{r},t)=E(\mathbf{r},t)$ is electric component and $\nabla\times A(\mathbf{r},t)=B(\mathbf{r},t)$ is the magnetic component of the optical field. The operators $a_{k}$ and $a_{k}^{\dagger}$ are known as photon annihilation and creation operators for the $k$th mode and obey the canonical commutation relations $[a_{k},a_{l}^{\dagger}]=\delta_{kl}$. In the presence of sources these operators become slowly varying functions of time or one of spatial coordinates.

In the following we consider a finite number $n$ of optical modes. We introduce a column vector of field operators $\mathbf{a}=(a_{1},...a_{n},a_{1}^{\dagger},...,a_{n}^{\dagger})^{T}$, and a column vector of modal functions $\mathbf{f}(\mathbf{r},t)=(f_{1}(\mathbf{r},t),...,f_{n}(\mathbf{r},t),f_{1}^{*}(\mathbf{r},t),...f_{n}^{*}(\mathbf{r},t))^{T}$. Now Eq. (1) can be rewritten in a matrix form as $A(\mathbf{r},t)=\mathbf{f}^{\dagger}(\mathbf{r},t)\mathbf{a}$, and the commutation relations as $\mathbf{a}\mathbf{a}^{\dagger}-\left(\mathbf{a}^{\dagger T}\mathbf{a}^{T}\right)^{T}=\mathbf{K}$, where

[TABLE]

and $\mathbb{I}$ is the $n\times n$ unit matrix. Here and below we adopt the bold font for $2n\times 2n$ matrices and $2n\times 1$ column vectors, the uppercase letters being reserved for matrices and the lowercase ones for vectors. The “blackboard bold” font is used for $n\times n$ matrices, while $n\times 1$ column vectors are marked by a bar above a lowercase letter. We introduce here a “complex symplectic matrix formalism”, where a Hermitian conjugation applied to a matrix of operators means a transposition of the matrix and a Hermitian conjugation of its elements. Transposition applied to such a matrix means a transposition of the matrix without any effect to the elements. This formalism differs from the traditional real symplectic formalism Simon et al. (1994); Ferraro et al. (2005); Adesso et al. (2014); Weedbrook et al. (2012) in two aspects. First, our variables are non-Hermitian operators of photon creation and ahhihilation and not positions and momenta. Such an approach is mentionned in Refs. Simon et al. (1994); Ferraro et al. (2005); Adesso et al. (2014), though not used for practical calculations. Second, all relations are expressed in a compact matrix form, simplifying the calculation, while in the traditional approach many important relations are expressed via matrix elements.

II.2 General Gaussian unitary transformation

Interaction between the modes in a nonlinear optical process results in a transformation of the field operators. The case of a linear in $\mathbf{a}$ transformation is called a Gaussian unitary transformation and has the form

[TABLE]

where $\mathbf{S}$ is a complex matrix. The conservation of the commutation relations requires $\mathbf{S}\mathbf{K}\mathbf{S}^{\dagger}=\mathbf{K}$, i.e. that the matrix $\mathbf{S}$ is complex symplectic, and has the following structure Simon et al. (1994); Ferraro et al. (2005); Adesso et al. (2014):

[TABLE]

where the complex $n\times n$ matrices $\mathbb{S}_{0}$ and $\mathbb{S}_{I}$, known as Bogoliubov transfromation matrices, satisfy the relation $\mathbb{S}_{0}\mathbb{S}_{0}^{\dagger}-\mathbb{S}_{I}\mathbb{S}_{I}^{\dagger}=\mathbb{I}$ and the matrix $\mathbb{S}_{0}\mathbb{S}_{I}^{T}$ is complex symmetric.

For each symplectic matrix $\mathbf{S}$ there exists a unitary operator $\mathcal{U}$, defined up to phase, producing the corresponding transformation of the field operators Simon et al. (1994)

[TABLE]

For two subsequent transformations we have $\mathbf{S}=\mathbf{S}_{2}\mathbf{S}_{1}$ and $\mathcal{U}=\mathcal{U}_{2}\,\mathcal{U}_{1}$.

Any unitary operator can be written as $\mathcal{U}=\exp(-i\mathcal{H})$, where the Hermitian operator $\mathcal{H}$ is called transformation generator. For a Gaussian unitary transformation its generator is a polynomial of $\mathbf{a}$ of the second order at maximum. For the sake of simplicity we omit the linear part of this dependence, which is rather trivial, and consider the generators being quadratic forms in the field operators: $\mathcal{H}=\frac{1}{2}\mathbf{a}^{\dagger}\mathbf{H}\mathbf{a}$, where the Hermitian matrix $\mathbf{H}$, which we call “transformation generator matrix”, has the structure

[TABLE]

with $\mathbb{H}_{0}$ Hermitian and $\mathbb{H}_{I}$ complex symmetric.

The symplectic matrix, Eq. (4), can also be expressed via its generator Adesso et al. (2014), which is simply related to the matrix of the quadratic form: $\mathbf{S}=\exp(-i\mathbf{K}\mathbf{H})$. Thus, any symplectic matrix can be written as

[TABLE]

The dimension of the group, created by symplectic matrices of the form Eq. (4) or Eq. (7) is $2n^{2}+n$. This is also the number of linearly independent generators of the form $\mathbf{K}\mathbf{H}$. We note that any pair of $n\times n$ matrices, a Hermitian $\mathbb{H}_{0}$ and a complex symmetric $\mathbb{H}_{I}$, define an $n$-dimensional Gaussian unitary transformation. Thus, parametrisation of a Gaussian unitary transformation in terms of these matrices is much simpler than in terms of the Bogoliubov transformation matrices $\mathbb{S}_{0}$ and $\mathbb{S}_{I}$, which should satisfy additional equations involving both matrices, as indicated above.

II.3 Passive Gaussian unitary transformation

An important subclass of Gaussian unitary transformations contains the transformations, preserving the total number of photons $\mathbf{a}^{\dagger}\mathbf{a}$. Such transformations correspond physically to mixing different modes on a multiport interferometer and are known as passive Gaussian unitary transformations Simon et al. (1994); Ferraro et al. (2005); Weedbrook et al. (2012); Adesso et al. (2014). The passive transformations necessarily have $\mathbb{H}_{I}=\mathbb{S}_{I}=0$, the matrix $\mathbb{S}_{0}$ becoming unitary: $\mathbb{S}_{0}=\exp(-i\mathbb{H}_{0})$. The symplectic matrix

[TABLE]

is also unitary in this case. The corresponding transformation generator can be obtained by writing

[TABLE]

where $\bar{a}=(a_{1},...a_{n})^{T}$ is the column vector of the annihilation operators of the modes. The corresponding unitary operator in the state space is (up to phase)

[TABLE]

The dimension of this (compact) subgroup of transformations is $n^{2}$.

II.4 Mode-wise squeezing

Another important subclass of Gaussian unitary transformations consists of single-mode squeezing of all modes with the squeezing parameters $\{r_{1},...,r_{n}\}$, some of which may be zero Simon et al. (1994); Ferraro et al. (2005); Weedbrook et al. (2012); Adesso et al. (2014). For definiteness we accept that all squeezing parameters are non-negative and sorted in the decreasing order. They can be written in a form of a positive diagonal matrix $\mathbb{R}=\mathrm{diag}\{r_{1},...,r_{n}\}$. We define the mode-wise squeezing transformation as a Gaussian unitary transformation with $\mathbb{H}_{0}=0$ and $\mathbb{H}_{I}=i\mathbb{R}$. The symplectic matrix for this class of transformations is

[TABLE]

The corresponding Hamiltonian can be obtained by writing

[TABLE]

The corresponding unitary operator in the state space is

[TABLE]

The dimension of this (noncompact) subgroup of transformations is $n$.

II.5 Bloch-Messiah reduction

The central mathematical procedure for determining the modes of squeezing of a multimode optical field is the decomposition of an arbitrary Gaussian unitary transformation into two passive transformations and one mode-wise squeezing transformation:

[TABLE]

where the matrices $\mathbb{V}$ and $\mathbb{Q}$ are unitary and the matrix $\mathbb{R}$ is positive diagonal. This decomposition is known as Bloch-Messiah reduction. It was introduced by Bloch and Messiah for fermions Bloch and Messiah (1962) and later generalized to bosons Balian et al. (1965). The physical meaning of this procedure is the possibility of realizing any Gaussian unitary transformation by means of two multiport interferometers and a number of single-mode squeezers Braunstein (2005). It should be noted that in practice a passive Gaussian transformation can be realized not by a multiport interferometer, but by a change of modal basis, corresponding to a proper shaping of the local oscillator used in the homodyne measurement of the field.

Decomposition defined by Eq. (24) can be written in the following form, used by Braunstein Braunstein (2005):

[TABLE]

which means a simultaneous SVD of two Bogoliubov transformation matrices.

The corresponding expression for the unitary operator in the state space is

[TABLE]

where the Hermitian matrices $\mathbb{H}_{V}$ and $\mathbb{H}_{Q}$ are the generators of the corresponding unitary matrices in Eq. (24): $\mathbb{V}=e^{-i\mathbb{H}_{V}}$, $\mathbb{Q}=e^{-i\mathbb{H}_{Q}}$.

II.6 Gaussian states

A fundamental property of a Gaussian transformation is that it transforms a Gaussian state into another Gaussian one. Any $n$-mode state with a density operator $\rho$ is characterized by its Wigner function

[TABLE]

where $\bar{\alpha}=(\alpha_{1},...,\alpha_{n})^{T}$ and $\bar{\lambda}=(\lambda_{1},...,\lambda_{n})^{T}$ are two column-vectors of complex variables and $\boldsymbol{\alpha}=(\bar{\alpha}^{T},\bar{\alpha}^{\dagger})^{T}$.

A state is called Gaussian if its Wigner function is a multivariate Gaussian distribution Simon et al. (1994); Ferraro et al. (2005); Adesso et al. (2014):

[TABLE]

where $\boldsymbol{\alpha}_{0}=\langle\mathbf{a}\rangle$ is the column vector of mean field and

[TABLE]

is the $2n\times 2n$ complex covariance matrix.

If a Gaussian state with mean $\boldsymbol{\alpha}_{0}$ and covariance matrix $\boldsymbol{\Sigma}$ undergoes a Gaussian transformation with matrix $\mathbf{S}$, then the state at the output is a Gaussian state with the mean $\mathbf{S}\boldsymbol{\alpha}_{0}$ and the covariance matrix $\boldsymbol{\Sigma}^{\prime}=\mathbf{S}\boldsymbol{\Sigma}\mathbf{S}^{\dagger}$ Simon et al. (1994).

The complex covariance matrix $\boldsymbol{\Sigma}$ has the following structure

[TABLE]

where the Hermitian matrix

[TABLE]

is the phase-insensitive covariance (also known as coherency matrix), while the complex symmetric matrix

[TABLE]

is the phase-sensitive covariance of the field (also known as anomalous correlator Kilin (1989)).

In the next section we will see how the developed formalism leads to a natural determination of squeezing eigenmodes in the process of PDC.

III Squeezing eigenmodes

III.1 Determining the squeezing eigenmodes by Bloch-Messiah reduction

In many nonlinear optical experiments, such as PDC and four-wave mixing (FWM), a squeezed light is generated by amplification of vacuum fluctuations. The formalism of Bloch-Messiah reduction, when applied to the vacuum field at the input, gives a possibility to define a set of modes at the output, such that each mode of the field is in a squeezed vacuum state with the squeezed quadrature along the same direction in the phase space, or vacuum. Indeed, the field transformation Eq. (3) can be written with the help of Eq. (24) as

[TABLE]

where $\mathbf{a}_{vac}=\mathbf{Q}^{\dagger}\mathbf{a}$ and the unitary matrices $\mathbf{V}$ and $\mathbf{Q}^{\dagger}$ are the leftmost and the rightmost ones in the right hand side of Eq. (24), defined as

[TABLE]

The modes characterized by operators $\mathbf{b}=\mathbf{V}^{\dagger}\mathbf{a}^{\prime}$ with the modal functions $\mathbf{g}(\mathbf{r},t)=\mathbf{V}^{\dagger}\mathbf{f}(\mathbf{r},t)$ are generally known as “modes of squeezing” or “squeezed modes”. We shall call them “squeezing eigenmodes” and the corresponding diagonal values of the matrix $\mathbb{R}$ “squeezing eigenvalues”. These modes are important for considering encoding quantum information into continuous variables of the optical field. For example, in homodyne detection, for observing the maximal squeezing, the modal function of the local oscillator should match the mode with the maximal squeezing eigenvalue $r_{1}$.

III.2 Determining the squeezing eigenmodes by squeezing matrix

An alternative way of determining the squeezing eigenmodes can be obtained from Eq. (26), which we rewrite as Ma and Rhodes (1990)

[TABLE]

where $\Gamma=\mathbb{V}\mathbb{R}\mathbb{V}^{T}$ is a compex symmetric matrix, called “the squeezing matrix” and we have used the identity

[TABLE]

following from the properties of a passive Gaussian transformation, described in Sec. II.3. When the operator $\mathcal{U}$ acts on a vacuum field, as in the case of unseeded PDC, two rightmost factors in the right hand side of Eq. (35) have no effect, and the resulting state is determined only by the third factor, dependent on $\Gamma$.

Bennink and Boyd Bennink and Boyd (2002) considered the Takagi factorization of the squeezing matrix $\Gamma=\mathbb{VRV}^{T}$ and called the modes defined by the columns of $\mathbb{V}$ the “eigenmodes of the squeezing”. We see from the above that the Takagi factorization of $\Gamma$ results in the same matrix $\mathbb{V}$ as the Bloch-Messiah reduction, and therefore this way of determining the squeezing eigenmodes is equivalent to that of Sec. III.1.

III.3 Determining the squeezing eigenmodes by covariance matrix

One more alternative way for determining the squeezing eigenmodes is based on the diagonalization of the covariance matrix. The state at the input of the nonlinear crystal in the case of unseeded PDC is vacuum, which is a Gaussian state with zero mean and the covariance matrix $\frac{1}{2}\mathbf{I}$. In accord with Sec. II.6 the output field has zero mean and the covariance matrix

[TABLE]

where

[TABLE]

The blocks of the covariance matrix, defined by Eq. (30) are $\Sigma^{\prime}_{0}=\frac{1}{2}\mathbb{V}\cosh(2\mathbb{R})\mathbb{V}^{\dagger}$ and $\Sigma^{\prime}_{I}=\frac{1}{2}\mathbb{V}\sinh(2\mathbb{R})\mathbb{V}^{T}$. Note, that the matrix $\mathbb{V}$ obtained from the Takagi factorization of $\Sigma^{\prime}_{I}$ fits also the spectral decomposition of $\Sigma^{\prime}_{0}$, but the inverse is not true in general.

We see thus that the squeezing eigenmodes can be obtained from the Takagi factorization of the phase-sensitive covariance $\Sigma^{\prime}_{I}$ of the output field, similar to the approach of Shapiro and Shakeel Shapiro and Shakeel (1997). In the latter approach the squeezing eigenmodes are obtained from a diagonalization of a real symmetric covariance matrix, including both the phase-insensitive and the phase-sensitive covariances. As we have shown above the same result can be achieved by a diagonalization of the phase-sensitive covariance alone. However this matrix is complex symmetric, not real symmetric as in Ref. Shapiro and Shakeel (1997).

As we have seen from the two last sections, the matrix of squeezing eigenmodes $\mathbb{V}$ can be obtained by the Takagi factorization of either the squeezing matrix $\Gamma$ determined by the physical model of PDC, or the phase-sensitive covariance matrix $\Sigma^{\prime}$, directly measurable by the homodyne technique. Takagi factorization can be rather easily calculated numerically when all the squeezing eigenvalues are different. This case is typical for degenerate OPA Wasilewski et al. (2006); Lvovsky et al. (2007) and optical parametric oscillators Patera et al. (2010); Jiang et al. (2012); Arzani et al. (2018). As we will see in Sec. IV, in the case of twin beams each squeezing eigenvalue has a multiplicity of at least two and the squeezing eigenmodes are defined with some degree of freedom. Takagi factorization in this case requires computing a balancing matrix and its square root Cariolaro and Pierobon (2016b), which introduces a higher level of complexity, especially in a multi-dimensional case. We show below how this complexity may be reduced in many practically important cases.

III.4 Frequency eigenmodes

To show how the formalism developed above is related to generation of squeezed light in PDC, we limit our consideration to the case, where the generated photons differ in one variable only, the frequency, i.e. to the case of type-I collinear PDC. However, all obtained results will be valid when the photons differ in any combination of frequency, direction and polarization with a proper relabeling of the modes.

We consider PDC in a nonlinear $\chi^{(2)}$ crystal with a pulsed plane-wave pump of central frequency $\omega_{p}$. We assume that the pump wave is strong enough and is undepleted. A coordinate system is chosen with the $z$-axis in the direction of the pump wave propagation and with the origin at the input edge of the crystal, see Fig. 1. The pump is considered as a classical wave with the (time-dependent) amplitude $E_{p}$, the wave vector $k_{p}$, and the frequency $\omega_{p}$. The downconverted wave may have a broad spectrum of frequencies $\omega=\omega_{0}+\Omega$ around the central frequency $\omega_{0}=\omega_{p}/2$, with the corresponding wave vector $k_{0}$.

The down-converted wave is described by the positive-frequency operator $E^{(+)}(t,z)$ normalized to photon-flux units, which can be decomposed into Fourier components as

[TABLE]

where $a(\Omega,z)$ is the photon annihilation operator with the frequency $\omega_{0}+\Omega$ and the longitudinal coordinate $z$.

Another operator, $\epsilon(\Omega,z)$, is defined by the relation Kolobov (1999)

[TABLE]

where $k(\Omega)$ is the wave vector of the down-converted light in the crystal, corresponding to the frequency $\omega_{0}+\Omega$. Operator $\epsilon(\Omega,z)$ is convenient for the description of the nonlinear interaction inside the crystal and is a quantum-mechanical analog of the classical slowly-varying amplitude Boyd (2008). The point $z_{0}$, where $a(\Omega,z_{0})=\epsilon(\Omega,z_{0})$, can be placed anywhere inside the crystal. Variation of this point results in a phase factor for the slowly varying amplitude. The transformation Eq. (40) can be alternatively considered as passage to the interaction picture Wasilewski et al. (2006); Lvovsky et al. (2007). This passage is dependent on the reference frame, and the point $z_{0}$ can be viewed as the origin of the “interaction picture reference frame”. We recall that the “working reference frame” has its origin at the crystal input, so that the field evolution in the crystal is considered between the points $z=0$ and $z=L$, where $L$ is the crystal length. Such a choice of the working reference frame simplifies the formulas. We analyze below several choices for $z_{0}$, including $z_{0}=0$ and $z_{0}=L/2$, used in the literature.

The classical pump field is decomposed into Fourier components as

[TABLE]

where $k_{p}(\Omega)$ is the wave vector of the pump at the frequency $\omega_{p}+\Omega$ and $\mathcal{E}^{\text{in}}(\Omega)$ is the Fourier transform of the pump wave at the crystal input at $z=0$. We note that for the type-I phase-matching the polarization of the pump wave is orthogonal to that of the generated light and the refractive index determining the dependence $k_{p}(\Omega)$ is different from that determining the dependence $k(\Omega)$.

Instead of working with continuous functions of frequency, we prefer to define discrete modes at frequencies, equidistantly spaced by the value $2\pi/T$, which is determined by the duration $T$ of the time interval where the quantum field is considered. We limit our consideration to the frequency band $[\omega_{0}-\Delta\Omega,\omega_{0}+\Delta\Omega]$, where the photons are mainly generated, and call the photons with positive (negative) detunings signal (idler) ones. We split each band into $m$ modes, and write the discrete frequencies as

[TABLE]

where $l$ runs from 1 to $n=2m$.

For these discrete modes we have a set of annihilation operators $a_{l}(z)=a(\Omega_{l},z)$ and a set of slowly-varying amplitudes $\epsilon_{l}(z)=\epsilon(\Omega_{l},z)$. Composing the vectors of length $2n$ we obtain $\mathbf{a}(z)=(a_{1}(z),...a_{n}(z),a_{1}^{\dagger}(z),...,a_{n}^{\dagger}(z))^{T}$, and $\boldsymbol{\xi}(z)=(\epsilon_{1}(z),...\epsilon_{n}(z),\epsilon_{1}^{\dagger}(z),...,\epsilon_{n}^{\dagger}(z))^{T}$.

The evolution of the downconverted light in the crystal is described by the following equation Wasilewski et al. (2006); Bennink and Boyd (2002):

[TABLE]

where $\mathcal{E}_{j+l}^{\text{in}}=\mathcal{E}^{\text{in}}(\Omega_{j}+\Omega_{l})$, and the first term in the right-hand side corresponds to dispersive propagation, while the second term corresponds to nonlinear interaction with the coupling constant $\sigma$, proportional to the nonlinear susceptibility of the crystal. In this nonlinear interaction a photon of the pump wave with the frequency $\Omega_{j}+\Omega_{l}$ is annihilated and converted into two photons with frequencies $\Omega_{j}$ and $\Omega_{l}$ respectively.

For slowly-varying amplitudes we obtain from Eqs. (40), (III.4)

[TABLE]

where

[TABLE]

is the Fourier component of the pump field at $z=z_{0}$, and $\Delta_{jl}=k_{p}(\Omega_{j}+\Omega_{l})-k(\Omega_{j})-k(\Omega_{l})$ is the phase mismatch of the corresponding modes. Equation (44) can be rewritten in a compact matrix form as

[TABLE]

where the matrix $\mathbf{F}(z)$ is given by

[TABLE]

with the complex symmetric matrix $\mathbb{F}(z)$ defined as

[TABLE]

Equation (46) can be written in a Hamiltonian form $\partial_{z}\boldsymbol{\xi}(z)=-i[\boldsymbol{\xi}(z),\mathcal{F}(z)]$ with a $z$-dependent Hamiltonian

[TABLE]

where $\bar{\epsilon}=(\epsilon_{1},...\epsilon_{n})^{T}$.

Solution of Eq. (46) can be written in the form of a $T$-exponent either for a symplectic matrix, transforming $\boldsymbol{\xi}(0)$ to $\boldsymbol{\xi}(L)$:

[TABLE]

or for a unitary operator

[TABLE]

where $\mathcal{T}$ denotes the $z$-ordering superoperator, placing the operators (matrices) with higher $z$ to the left in the Taylor series of the exponential. Comparing Eq. (51) with the definition of the transformation generator in Sec. II.2, we conclude that the latter is given by $\mathcal{H}=\mathcal{F}L$ in the case of $z$-independent Hamiltonian, which can be met, e.g. in the case of perfect phasematching for all modes, where $\Delta_{jl}=0$. This fact explains the widely used referring to the transformation generator $\mathcal{H}$ as “Hamiltonian”. However, in the general case the Hamiltonian of the field is $z$-dependent and is defined by Eq. (49).

III.5 Magnus expansion

Both $T$-exponents, Eqs. (50) and (51), can be represented as Magnus expansions Magnus (1954); Blanes et al. (2009)

[TABLE]

where $\boldsymbol{\Omega}_{k}$ is a $2n\times 2n$ matrix proportional to $|\sigma|^{k}$. The first term in Eq. (52) is

[TABLE]

and the higher order terms are expressed via integrals over nested commutators of $\mathbf{F}(z)$ with itself at different spatial points. Equivalently, $\boldsymbol{\Xi}_{k}$ is an operator proportional to $|\sigma|^{k}$, and the first term in Eq. (53) is

[TABLE]

while the higher order terms are expressed via integrals over nested commutators of $\mathcal{F}(z)$ with itself at different spatial points.

First order Magnus approximation to the symplectic matrix is denoted as $\mathbf{S}_{1}$ and obtained by leaving the first term only in the exponent of Eq. (52), while a similar approximation to the evolution operator is denoted $\mathcal{U}_{1}$ and obtained by leaving the first term only in the exponent of Eq. (53). Let us show that $\mathcal{U}_{1}$ is a Gaussian transformation with the corresponding symplectic matrix $\mathbf{S}_{1}$. Substituting Eq. (49) into Eq. (55) we obtain

[TABLE]

corresponding to a Gaussian transformation with the symplectic matrix $e^{\mathbf{K}^{2}\boldsymbol{\Omega}_{1}}=e^{\boldsymbol{\Omega}_{1}}$, which is exactly $\mathbf{S}_{1}$.

The symplectic matrix $\mathbf{S}_{1}$ is determined by its generator, which we denote by $\mathbf{H}^{[1]}$, as $\mathbf{S}_{1}=\exp\left(-i\mathbf{K}\mathbf{H}^{[1]}\right)$. The generator is determined according to Eq. (6) by two $n\times n$ matrices: $\mathbb{H}_{0}^{[1]}$, which is Hermitian, and $\mathbb{H}_{I}^{[1]}$, which is complex symmetric, as

[TABLE]

We have from Eqs. (54), (47) $\mathbb{H}_{0}^{[1]}=0$, and

[TABLE]

wherefrom with the help of Eq. (48) we obtain

[TABLE]

From Eq. (56) we have

[TABLE]

Comparing this expression with Eq. (35), we arrive at the conclusion that in the first order of Magnus expansion the squeezing matrix is $\Gamma=-i\mathbb{H}_{I}^{[1]}$. As consequence, the squeezing eigenmodes are determined by the Takagi factorization $-i\mathbb{H}_{I}^{[1]}=\mathbb{VRV}^{T}$, where $\mathbb{V}$ is unitary. We note that the Takagi factorization for $i\mathbb{H}_{I}^{[1]*}$ reads $i\mathbb{H}_{I}^{[1]*}=\mathbb{V}^{*}\mathbb{RV}^{\dagger}$. When the matrices $\mathbb{V}$ and $\mathbb{R}$ are found, the Bloch-Messiah reduction in the first order of Magnus expansion can be easily written in the form

[TABLE]

which is a special case of Eq. (24) with $\mathbb{Q}=\mathbb{V}$. Note that in this special case the coincidence of the input and the output eigenmodes reflects certain symmetry of the first order of Magnus expansion which does not hold in the higher orders.

The unitary evolution operator in the state space is given by Eq. (26) with $\mathbb{H}_{Q}=\mathbb{H}_{V}$. When this operator acts on the vacuum, it produces a multimode squeezed state

[TABLE]

In the new basis, defined as

[TABLE]

the transformation of the field operators has a form of mode-wise squeezing:

[TABLE]

The new modal functions are $\mathbf{f}^{\prime}(\mathbf{r},t)=\mathbb{V}^{\dagger}\mathbf{f}(\mathbf{r},t)$. Note, that the field operator is invariant with respect to the choice of the modal basis: $A(\mathbf{r},t)=\mathbf{f}^{\prime}{}^{\dagger}(\mathbf{r},t)\boldsymbol{\xi}^{\prime}=\mathbf{f}^{\dagger}(\mathbf{r},t)\boldsymbol{\xi}$.

In the present work we limit our consideration to the first order of the Magnus expansion, which has been found a good approximation for not very high squeezing, below 12 dB Christ et al. (2013); Lipfert et al. (2018). An analytic treatment of higher orders of the Magnus expansion in PDC can be found in Refs. Quesada and Sipe (2014, 2015); Lipfert et al. (2018).

III.6 Reduction of the Takagi factorization to a real symmetric spectral decomposition

As we have seen above, in the first order of the Magnus expansion the squeezing eigenmodes are defined by the Takagi factorization applied to the complex symmetric squeezing matrix $\Gamma=-i\mathbb{H}_{I}^{[1]}$, whose matrix elements are determined by Eq. (59). Takagi factorization can always be realized for a complex symmetric matrix by performing a singular value decomposition and then computing a balancing matrix Cariolaro and Pierobon (2016a, b). However, this procedure can be significantly simplified under certain conditions. We notice that the elements of the squeezing matrix can be made real if two conditions are satisfied: (i) the origin of the interaction picture frame is chosen at the center of the crystal, $z_{0}=L/2$, and (ii) the pump pulse at this point is transform-limited, i.e. all elements of the vector $\mathcal{E}_{j+l}^{0}$ have the same phase $\varphi_{0}$. Then, as is easily seen from Eq. (59), all elements of the squeezing matrix have the phase $\varphi=\varphi_{0}+\arg(\sigma)-\pi/2$. This phase can be removed by a trivial transformation $\bar{\epsilon}\to\bar{\epsilon}e^{i\varphi/2}$ in Eq. (46) which makes the squeezing matrix real.

When the squeezing matrix is real symmetric it can be represented in the form of spectral decomposition

[TABLE]

where $\mathbb{O}$ is a real orthogonal matrix and $\Lambda$ is the real diagonal matrix of eigenvalues. If all eigenvalues of $\Gamma$ are non-negative, i.e. the squeezing matrix is positive semidefinite Horn and Johnson (1985), Eq. (76) is a Takagi factorization, the columns of $\mathbb{O}$ are the modal functions of the modes of squeezing, and the eigenvalues of $\Gamma$ are the squeezing parameters. In the opposite case, when some of the eigenvalues of $\Gamma$ are negative, the Takagi factorization can be reconstructed from the spectral decomposition in the following way. Let $\bar{e}_{j}=(0,...,1,...0)^{T}$ be the $1\times n$ column vector with the $j$th element equal to 1 and all others zero. Then $\mathbb{O}\bar{e}_{j}=O_{j}$ is the eigenvector corresponding to the $j$th eigenvalue $\lambda_{j}$. We define a set of column vectors $V_{j}$ as follows: $V_{j}=O_{j}$ if $\lambda_{j}\geq 0$ and $V_{j}=iO_{j}$ if $\lambda_{j}<0$. We build a square matrix $\mathbb{V}=(V_{1},...,V_{n})$ from the columns $V_{j}$. It is easy to verify that this matrix is unitary, since $\mathbb{O}$ is orthogonal. We define also $r_{j}=|\lambda_{j}|$ and build a diagonal matrix $\mathbb{R}$ with $r_{j}$ at diagonal. Now we can rewrite Eq. (76) as

[TABLE]

which is a Takagi factorization for $\Gamma$. Thus, in the case of real squeezing matrix the modal functions of the squeezing eigenmodes, defined by the columns of $\mathbb{V}$, are either purely real or purely imaginary, which greatly simplifies their analysis and visualisation.

The price to be paid for the possibility to work with real modal functions is, as shown above, the necessity to have the pump pulse transform-limited at the center of the crystal. In a physical experiment the pump pulse is typically generated as transform-limited by the laser, and is therefore transform-limited at the crystal input. During its propagation inside the crystal the pulse experiences dispersion determined by the dependence $k_{p}(\Omega)$ in Eq. (45). As a result, at the point $z_{0}=L/2$ it acquires the phase $k_{p}(\Omega)L/2$. In the quadratic dispersion approximation we can write

[TABLE]

where $k_{p0}$, $k_{p0}^{\prime}$ and $k_{p0}^{\prime\prime}$ are the value and the two derivatives of $k_{p}(\Omega)$ at $\Omega=0$. Thus, the phase of the pump pulse at the crystal center is composed of three major terms: (i) constant phase shift due to optical oscillations at the carrier frequency, which can be included into the phase $\varphi_{0}$, (ii) absolute group delay of the pump pulse, which can be put to zero by adjusting the time origin, and (iii) the frequency chirp, which is equivalent to pulse spread in the temporal domain. The latter phase cannot be removed by a simple mathematical transformation and its removal has physical meaning. First, it can be neglected if the pump bandwidth $\Omega_{p}$ is small enough so that $\Omega_{p}\ll 2\sqrt{2}(k_{p0}^{\prime\prime}L)^{-1/2}$. Second, it can be compensated by pre-chirping the pump pulse with the frequency chirp $-k_{p0}^{\prime\prime}L\Omega^{2}/4$. Since most nonlinear crystals have positive dispersion ($k_{p0}^{\prime\prime}>0$), this pre-chirp should be negative and thus require an active phase control.

In the examples of the squeezing eigenmodes presented below in Sec. V we suppose that one of these conditions holds and the modal functions of squeezing eigenmodes, defined in terms of the slowly-varying amplitudes $\epsilon(\Omega)$, are either real or imaginary. For the choice of $z_{0}=L/2$ the slowly-varying amplitude coincides with the full field $a(\Omega)$ at the crystal center, and at the crystal output the latter operator has additional phase determined by Eq. (40). This phase includes the three terms shown in Eq. (78) for the pump: constant phase shift, group delay and chirp. The full analysis of the output modal functions is outside of the scope of the present treatment, which is limited to the squeezing eigenvalues and eigenmodes at the crystal center. We note that recently the eigenmodes of a similar parametric nonlinear process – sum-frequency generation – have been found for the case of pre-chirped pump pulse Patera et al. (2018).

III.7 Summary

In this section we have shown how the squeezing eigenmodes and the corresponding squeezing eigenvalues can be obtained from the Takagi factorization of the squeezing matrix of the multimode optical field. Further, we have shown how the squeezing matrix can be obtained from the physical model of PDC, and discussed a special case where this matrix is real. In the next section we apply this general theory to twin beams of light, having remarkable additional symmetries.

IV Twin beams of light

IV.1 Parametric downconversion in the case of twin beams

Up to this point our consideration of PDC related any case of phase-matching: degenerate and nondegenerate. In the following sections we apply the developed formalism to the specific case of twin beams and demonstrate some important regularities inherent in this case. Twin beams are generated in non-degenerate PDC where signal and idler photons are discriminated by frequency, direction, polarization or a combination of these variables. For the description of twin beams we need to consider a set of $m$ signal modes and $m$ idler modes which should be united in the formalism of the previous sections with the total number of $n=2m$ modes. The fact that in every photon pair the photons are always created in different sets of modes imposes additional constraints on the matrices $\mathbb{H}_{0}$ and $\mathbb{H}_{I}$, defining the Gaussian unitary transformation.

To be specific, we again limit our consideration to the case of twin beams discriminated by frequency, with obvious extensions to other cases. Twin beams generation is characterized by phase-matching conditions, providing emission of signal photons in a frequency band $[\omega_{s}-\Delta\Omega_{0},\omega_{s}+\Delta\Omega_{0}]$ and idler photons in a frequency band $[\omega_{i}-\Delta\Omega_{0},\omega_{i}+\Delta\Omega_{0}]$ such that the frequency gap between the highest idler frequency $\omega_{i}+\Delta\Omega_{0}$ and lowest signal frequency $\omega_{s}-\Delta\Omega_{0}$ is much higher than the width of the pump spectrum. In this case emission of two photons into one band (signal or idler) is practically impossible, because the sum of their frequencies would be outside of the spectrum of the pump. As consequence, the elements of the matrix $\mathbb{H}_{I}^{[1]}$, defined by Eq. (59), are almost zero when two indices belong to the same group of modes (signal or idler). With a good degree of approximation we can put them to zero exactly and write the matrix $\mathbf{H}^{[1]}$ defining the generator of the Gaussian transformation in the first-order Magnus expansion as

[TABLE]

where $J$ is a complex $m\times m$ matrix, known as joint spectral amplitude (JSA) for two photons generated in an elementary nonlinear process. This matrix in general possesses no properties of symmetry.

As indicated above, our consideration here is limited to the first order Magnus approximation and for simplicity we omit the superscripts and the subscripts indicating the first order of Magnus expansion. We adopt the name of “twin beam Gaussian transformation” for a unitary transformation $\mathcal{U}=e^{-i\mathcal{H}}$ with a generator $\mathcal{H}=\frac{1}{2}\mathbf{a}^{\dagger}\mathbf{H}\mathbf{a}$, where the matrix $\mathbf{H}$ has the structure defined by Eq. (79).

The main regularity inherent in twin beams is stated by the following theorem.

Theorem 1**.**

All squeezing eigenvalues of a twin beam Gaussian transformation have multiplicities of at least two.

Proof.

Applying the SVD to the matrix $-iJ$, we write $-iJ=CR_{J}D^{\dagger}$, where $C$ and $D$ are unitary $m\times m$ matrices and $R_{J}$ is a diagonal $m\times m$ matrix with non-negative elements. Now we can write the $n\times n$ matrix $\mathbb{H}_{I}$, defined according to Eq. (6), as

[TABLE]

where

[TABLE]

Therefore the full $2n\times 2n$ generator matrix defined by Eq. (79) is

[TABLE]

and the corresponding symplectic matrix is

[TABLE]

which corresponds to a decomposition into two-mode squeezers. Note, that this decomposition is not the Bloch-Messiah reduction, since the matrix $\mathbb{\tilde{R}}$ is not diagonal.

The true Bloch-Messiah reduction can be obtained by an additional rotation. For each pair of conjugated modes the sub-matrix of $\mathbb{\tilde{R}}$ is proportional to the Pauli matrix $\sigma_{x}$. The Takagi factorization for this matrix is $\sigma_{x}=xx^{T}$, where

[TABLE]

Defining the matrix $\mathbb{X}$ as a direct sum of $m$ matrices $x$ for each pair of conjugated modes, we rewrite Eq. (80) as

[TABLE]

where the diagonal matrix

[TABLE]

is the matrix of squeezing eigenvalues. Now, defining

[TABLE]

we interpret Eq. (85) as the Takagi factorization of the squeezing matrix, which is equivalent to the Bloch-Messiah decomposition in the form Eq. (72). It is clear from Eq. (86) that each squeezing eigenvalue is found twice in the diagonal of $\mathbb{R}$, and therefore has multiplicity of at least two. ∎

The columns of matrices $C$ and $D$ are known as modal functions of the signal and idler Schmidt modes of a photon pair generated in the elementary process of photon-pair creation Law et al. (2000); Grice et al. (2001). It follows from the proof of Theorem 1, that the two squeezing eigenmodes corresponding to the same eigenvalue can be constructed from the Schmidt modes in a straightforward way. Let us denote the $j$th columns of the matrices $C$ and $D$ by $C_{j}$ and $D_{j}$ respectively, similar to how it was done in Sec. III.6 for $V_{j}$. Both $C_{j}$ and $D_{j}$ are $m\times 1$ matrices (column vectors) and their corresponding singular value is $r_{j}$. From Eqs. (87), (84) we obtain a $m\times 2$ matrix containing two squeezing eigenmodes corresponding to the same eigenvalue $r_{j}$

[TABLE]

Thus, one squeezing eigenmode can be constructed by concatenating the signal Schmidt mode $C_{j}$ and the conjugated idler Schmidt mode $D_{j}$. Another squeezing eigenmode can be constructed in a similar way from $iC_{j}$ and $iD_{j}$.

Let us study the uniqueness of the squeezing eigenmodes constructed in this way. For simplicity we suggest that each singular eigenvalue $r_{j}$ of $-iJ$ has multiplicity 1. Then by the Autonne’s uniqueness theorem Horn and Johnson (1985) the eigenvectors $C_{j}$ and $D_{j}$ are defined up to arbitrary phase $\varphi_{j}$, i.e. $C_{j}e^{i\varphi_{j}}$ and $D_{j}e^{i\varphi_{j}}$ are eigenvectors of an alternative SVD of $-iJ$. It means that the general form of Eq. (88) is

[TABLE]

It follows from the last expression that the squeezing eigenmodes corresponding to a given eigenvalue with multiplicity 2 are defined up to an orthogonal rotation. This rotation corresponds to shifting the phases of the signal and idler Schmidt modes by the same angle.

The multiplicity of the eigenvalues and the lack of definiteness of the eigenfunctions can be understood from considering just two modes with photon annihilation operators $a$ and $b$. The following identity holds Milburn (1984):

[TABLE]

where the real number $r$ is the squeezing parameter Kolobov (1999) and we have introduced two other modes with photon annihilation operators $c=e^{i\pi/4}(a-ib)/\sqrt{2}$ and $d=e^{i\pi/4}(b-ia)/\sqrt{2}$. The left hand side of Eq. (95) represents the two-mode squeeze operator Caves and Schumaker (1985), while the right hand side of this equation represents a product of two single-mode squeeze operators. Thus, a two-mode squeezed state can be represented by a unitary rotation in the modal space as direct product of two single-mode squeezed states with the same squeezing parameter. An important and less known feature of this representation consists in invariance of the right hand side of Eq. (95) with respect to orthogonal rotations in the modal space, $c\to c\cos\phi+d\sin\phi$, $d\to-c\sin\phi+d\cos\phi$, where $\phi$ is an arbitrary angle. Such a rotation is similar to that described by Eq. (91).

The practical importance of Theorem 1 is connected to its general applicability when the signal and the idler waves are well separated. Being a consequence of a symmetry condition, the multiplicity of the eigenvalues can be used as a priori information for fitting and reconstruction of the field state at the measurement stage.

IV.2 Reduction of the Takagi factorization to a spectral decomposition for twin beams

We have seen that the Takagi factorization of the squeezing matrix $\Gamma=-i\mathbb{H}_{I}$ can be obtained by an SVD in a two times smaller space of JSA. An alternative method consists in finding a spectral decomposition for the Hermitian matrix (cf. Theorem 7.3.3. of Ref. Horn and Johnson (1985))

[TABLE]

which we call “associated squeezing matrix”. It can be shown by direct substitution that the vectors

[TABLE]

are eigenvectors of $\Gamma_{a}$ with the eigenvalues $r_{j}$ and $-r_{j}$ respectively. Comparing Eq. (97) with Eq. (88) we arrive at a conclusion that the modal functions of the squeezing eigenmodes can be obtained by finding the eigenvectors of $\Gamma_{a}$ and then complex conjugating the idler part and shifting the phase of the vector with the negative eigenvalue by $\pi/2$.

This method is especially fruitful in the case of real matrix $\Gamma=-i\mathbb{H}_{I}$, where we have a situation already met in Sec. III.6. If we consider a pump pulse which is transform-limited at the center of the crystal, then the squeezing matrix is real (for an appropriately chosen pump phase), as follows from Eq. (59) with $z_{0}=L/2$:

[TABLE]

The associated squeezing matrix $\Gamma_{a}$ in this case coincides with $\Gamma$ and is also real symmetric. The eigenvectors $U_{j}^{(+)}$ and $U_{j}^{(-)}$ can be chosen real in this case Horn and Johnson (1985). Note that the eigenvectors of a complex spectral decomposition, Eq. (96) are complex vectors defined up to an arbitrary phase, while the eigenvectors of the real spectral decomposition, Eq. (76), are real vectors defined up to a sign. It means that by passing from a complex to a real matrix we find the eigenmodes in a more definite and reproducible way.

Indeed, choosing $U_{j}^{(+)}$ and $U_{j}^{(-)}$ real, we can omit the complex conjugation of the idler part of these vectors. It follows that $U_{j}^{(+)}$ and $iU_{j}^{(-)}$ are the modal functions of the squeezing eigenmodes corresponding to the squeezing eigenvalue $r_{j}$ with the multiplicity 2. Note that the two modal functions found by this algorithm are almost uniquely defined (up to sign), one of them being purely real, the other one being purely imaginary. Of course, orthogonal rotations in the space of these two functions give us other possible modal functions of the squeezing eigenmodes.

V Example of twin beams discriminated by frequency

V.1 Nondegenerate phase matching: numerical solution

As an example of the theory developed in the previous sections we consider twin beams generated by PDC in a $2$ mm long sample of BBO (beta barium borate, $\beta$-BaB2O4) crystal. For an angle between the crystal optical axis and the pump wave-vector of \theta_{0}=$$$ the type-I phase matching is satisfied for pumping at 397.5 nm, the signal and idler frequency being around 677 nm and 963 nm, respectively. We consider a Gaussian pump pulse, which is transform limited at the center of the crystal with the intensity full width at half maximum \tau_{p}=129$fs, corresponding to a spectral width of$1.8nm. The squeezing matrix is given in this case by Eq. ([98](#S4.E98)) and is shown in Fig. [2](#S5.F2) as function of detunings\Omega$and$\Omega^{\prime}$. The signal and the idler bands shown in Fig. 2 are well separated, which corresponds to twin beams generation.

We solve numerically the eigenvalue problem for real symmetric matrix $\Gamma=-i\mathbb{H}_{I}$, using the Sellmeier equation for the dispersion of the BBO crystal. The squeezing eigenvalues, given by the absolute values of the eigenvalues of $\Gamma$, $r_{k}=|\lambda_{k}|$, are shown in Fig. 3. We see that each squeezing eigenvalue has multiplicity 2, as stated by Theorem 1, and its logarithm is a linear function of the eigenmode pair number, so that the eigenvalues can be approximated as

[TABLE]

with $q_{num}=0.8903$ and $l$ running from 0 to $\infty$.

The corresponding modal functions of the squeezing eigenmodes are shown in Fig. 4. Modal functions corresponding to the positive eigenvalues of $\Gamma$ are real. Modal functions corresponding to its negative eigenvalues are purely imaginary and their imaginary part is plotted. The structure of the modal functions corresponds to that prescribed by Eq. (88): they are concatenations of Schmidt modal functions for the signal and the idler beams.

The effective number of squeezing eigenmodes can be estimated by a measure similar to the Schmidt number, widely used in the photon-pair generation regime Law and Eberly (2004); Horoshko et al. (2012); Gatti et al. (2012):

[TABLE]

where we have used the approximate expression for the squeezing eigenvalues, Eq. (99). We see that as $q_{num}$ approaches unity the effective number of squeezing eigenmodes tends to infinity.

V.2 Nondegenerate phase matching: analytic Gaussian modeling

As has been pointed out before, the problem of finding the squeezing eigenvalues and the squeezing eigenmodes can be approached by the SVD of the JSA matrix $J$, rather than by the Takagi factorization of the squeezing matrix $\Gamma$. In this section we follow this approach. The JSA matrix, defined by Eq. (79), is given by the top right $m\times m$ block of the matrix $\mathbb{H}_{I}^{[1]}$, determined by Eq. (59). We set $z_{0}=L/2$, as in the previous section, but here we consider a pump pulse transform-limited at the crystal input face, which is a typical experimental situation, so that at the crystal center it acquires a frequency-dependent phase determined by Eq. (45).

One advantage of the JSA approach consists in the economy of computational resources: the JSA matrix can be limited only to the spectral regions where the signal and the idler amplitudes are nonnegligible, while the squeezing matrix is mainly filled by zeros in the case of well-separated twin beams. Another advantage is the possibility to find the squeezing eigenvalues and the squeezing eigenmodes approximately by replacing the JSA with a double-Gaussian function and applying the well-known spectral decomposition of a real symmetric double-Gaussian kernel Grice et al. (2001); Wasilewski et al. (2006); Lvovsky et al. (2007); Patera et al. (2010); Horoshko et al. (2012)

[TABLE]

where $-1<q<1$, $p=\sqrt{1-q^{2}}$ and $h_{k}(x)=\left(2^{k}k!\sqrt{\pi}\right)^{-\frac{1}{2}}H_{k}(x)e^{-x^{2}/2}$ is the Hermite-Gauss function, $H_{k}(x)$ being the Hermite polynomial. Equation (101) is obtained by multiplying both sides of the Mehler’s formula for Hermite polynomials Mehler (1866) by $e^{-x^{2}/2-y^{2}/2}$. Since the Hermite-Gauss functions are orthonormal and complete on the Hilbert space, they represent eigenfunctions of the double-Gaussian kernel, while $pq^{k}$ are the corresponding eigenvalues.

Though Eq. (101) proved to be highly efficient for finding the eigenmodes and the eigenvalues of a degenerate PDC, in the nondegenerate case one needs an SVD for a complex double-Gaussian kernel, in general not symmetric. Such a decomposition can be obtained in the following form (see proof in Appendix):

[TABLE]

where $\mu,\nu,\eta$ and $\xi$ are real numbers, satisfying the relations $\mu,\nu>0$, $\sqrt{\mu\nu}>|\eta|$, necessary and sufficient for this kernel to be square-integrable. The parameters in the right-hand side of Eq. (V.2) are $\tau_{1}=\sqrt{uv/\nu}$, $\tau_{2}=\sqrt{uv/\mu}$, $\zeta=\eta\xi/(2uv)$, $p=\sqrt{1-q^{2}}$, and

[TABLE]

where $u=\sqrt{\mu\nu+\xi^{2}}$ and $v=\sqrt{\mu\nu-\eta^{2}}$. The phase $\theta_{k}$ is determined in Appendix. Note that the factor $v^{-\frac{1}{2}}$ in the right-hand side of Eq. (V.2) is the Hilbert-Schmidt norm of the kernel in its left-hand side. In Eq. (101) this norm is equal to 1. Another difference in the structure of two decompositions is the sign of the parameter $q$. In Eq. (101) this parameter can be negative, while in Eq. (V.2) it is always positive and is complemented by the phase $\theta_{k}$. In the limiting case $\xi=0$, $\mu=\nu=\sqrt{1+\eta^{2}}$, $\theta_{k}=\pi k(1-\eta/|\eta|)/2$, Eq. (V.2) reproduces Eq. (101).

Equation (V.2) shows that the left-singular functions of a complex double-Gaussian kernel are chirped Hermite-Gauss functions $\psi_{k}(x)=\sqrt{\tau_{1}}h_{k}(\tau_{1}x)e^{i\zeta\tau_{1}^{2}x^{2}+i\theta_{k}/2}$. It is straightforward to verify that these functions create a complete orthonormal set on the Hilbert space. The right-singular functions of this kernel are obtained from the left-singular ones by repacing $\tau_{1}$ by $\tau_{2}$.

In this section for the sake of analytic treatment we write the JSA function as a kernel depending on two continuous frequencies $\Omega_{1}$ and $\Omega_{2}$, but similar approach applies to the matrices obtained by passing to discrete frequency modes. To write the JSA in a double-Gaussian form, we make four following approximations.

The first approximation is the quadratic approximation for the dispersion law, which is a good approximation for a not too broadband PDC Caspani et al. (2010); Horoshko et al. (2012); Gatti et al. (2012). It has been applied to the pump in Eq. (78). Applying it to the downconverted light results in limiting the Taylor series of the phase mismatch function $\Delta(\Omega_{1},\Omega_{2})=k_{p}(\Omega_{1}+\Omega_{2})-k(\Omega_{1})-k(\Omega_{2})$ to the terms up to the second order in both frequencies:

[TABLE]

where $\Delta_{0}=k_{p0}-2k_{0}$ and we have introduced new variables $\Omega_{\pm}=\Omega_{1}\pm\Omega_{2}$. Here we denote by $k_{0}$, $k_{0}^{\prime}$ and $k_{0}^{\prime\prime}$ the value and the two derivatives of $k(\Omega)$ at $\Omega=0$, as we did for the pump in Sec. III.6.

The second approximation consists in disregarding the third term in the right hand side of Eq. (V.2) with respect to the second one, which can be done if $|\Omega_{+}|\ll 4|k_{p0}^{\prime}-k_{0}^{\prime}|/|2k_{p0}^{\prime\prime}-k_{0}^{\prime\prime}|$. For the considered example this means $|\Omega_{+}|\ll 1.04\omega_{0}$, while $|\Omega_{+}|$ is limited by the pump bandwidth $\Omega_{p}=2\sqrt{\log 2}/\tau_{p}\approx 0.0091\omega_{0}$, which means that this approximation is well justified. In this approximation the line in the frequency space corresponding to the perfect phase matching is represented by a parabola

[TABLE]

This parabola, corresponding to the middle of the phase-matched area in Fig. 2, intersects the $\Omega_{+}=0$ line at two points $\Omega_{-}=\pm 2\Omega_{s}$, where $\Omega_{s}=\sqrt{\Delta_{0}/k_{0}^{\prime\prime}}$ is the central detuning of the signal beam, $-\Omega_{s}$ being the central detuning of the idler one. In the vicinity of the point $\Omega_{+}=0$, $\Omega_{-}=2\Omega_{s}$ we introduce new variables $\delta\Omega_{+}=\Omega_{+}$, $\delta\Omega_{-}=\Omega_{-}-2\Omega_{s}$, corresponding to small deflections from the frequencies where the signal-idler coupling is maximal. Further, we disregard $\delta\Omega_{-}^{2}$ compared to $\Omega_{s}\delta\Omega_{-}$, which is the third approximation we make. This approximation is always valid for twin beams well-separated in frequency, when the bandwidth of one beam is much less than the frequency distance between the central frequencies of the beams.

In these three approximations the JSA function is obtained from Eqs. (59), (45) and (78) as

[TABLE]

where $g=\sigma LE_{0}e^{ik_{p0}L/2}$ is the coupling constant, $E_{0}$ being the peak pump amplitude, and we have introduced four characteristic times: the absolute group delay time of the pump $\tau_{pd}=k_{p0}^{\prime}L$, the pump spread time $\tau_{ps}=\sqrt{k_{p0}^{\prime\prime}L}$, the relative pump-signal group delay time $\tau_{d}=(k_{p0}^{\prime}-k_{0}^{\prime})L$, and $\tau_{s}=\sqrt{k_{0}^{\prime\prime}\Delta_{0}}L$. The latter time is proportional to the spread time of the signal pulse during its passage through the crystal $\sqrt{k_{0}^{\prime\prime}L}$ Horoshko et al. (2012); Gatti et al. (2012), but also depends on the phase-mismatch at degeneracy $\Delta_{0}$. It is straightforward to find that in the vicinity of the point $\Omega_{+}=0$, $\Omega_{-}=-2\Omega_{s}$ the signal-idler coupling is determined by the function $J(\Omega_{2},\Omega_{1})$, which is the transposition of the JSA defined by Eq. (106).

The fourth and the final approximation consists in making a replacement $\operatorname{sinc}(x)\approx e^{-x^{2}/(2\sigma_{0}^{2})}$, where $\sigma_{0}=1.61$ is chosen so that these two functions have the same width at half-maximum Grice et al. (2001); Wasilewski et al. (2006); Lvovsky et al. (2007); Patera et al. (2010); Horoshko et al. (2012). As a result, JSA in the variables $\delta\Omega_{1}=\Omega_{1}-\Omega_{s}$, $\delta\Omega_{2}=\Omega_{2}+\Omega_{s}$ takes the form $J(\Omega_{1},\Omega_{2})=g\tilde{J}(\Omega_{1},\Omega_{2})e^{i\tau_{pd}\left(\delta\Omega_{1}+\delta\Omega_{2}\right)/2+i\tau_{ps}^{2}\left(\delta\Omega_{1}^{2}+\delta\Omega_{2}^{2}\right)/4}$, where

[TABLE]

is the part of the kernel with a non-separable phase. This function is shown (in its absolute value) in Fig. 5, where it is compared to the exact one.

Comparing Eq. (V.2) to Eq. (V.2) and accepting for simplicity that $\arg(-ig)=0$, we find the modal functions of the signal and the idler Schmidt modes:

[TABLE]

where the parameters $\tau_{1}$ and $\tau_{2}$ are defined as in Eq. (V.2) with $\xi=\tau_{ps}^{2}/2$ and

[TABLE]

while the signal and the idler chirp rates are $\zeta_{1}=\xi(\nu+\eta)/(2uv)$ and $\zeta_{2}=\xi(\mu+\eta)/(2uv)$.

For the considered example $\tau_{1}=48$ fs, $\tau_{2}=60$ fs, $\zeta_{1}=0.0086$, $\zeta_{2}=-0.0063$, and $q=0.8681$. The modal function of the signal Schmidt mode $C_{1}(\Omega_{1})$ is shown in Fig. 6 together with the modal function obtained by a numerical SVD of JSA for two different crystal lengths. We see that the analytic solution follows closely the numerical one for a shorter crystal, with slight differences for a longer one.

We note also, that the variation of the phase of the modal function due to its chirp within the signal frequency band is much less than $\pi$ and therefore negligible for the considered example. This fact justifies the approach of the previous section, where the pump chirp was disregarded from the beginning.

It is interesting to analyze the influence of the crystal length $L$ on the shape of the Schmidt modes. As typical for PDC, growing crystal length results in a more restrictive phase matching and to narrower signal and idler bandwidths. This is clearly seen in Fig. 6. The absolute value of the chirp increases for a longer crystal, however, the maximal variation of the phase remains small because of decreased bandwidth. The dimensionless chirp parameters $\zeta_{1}$ and $\zeta_{2}$ allow us to estimate the importance of considering the complex squeezing matrix instead of a real one. In the considered example they are much less that 1. Asymptotically with growing $L$ they tend to length-independent constants, determined by the dispersive properties of the crystal material only, as can be seen from their definitions and Eq. (109).

We observe that the idler characteristic time $\tau_{2}$ is 25% higher than the signal characteristic time $\tau_{1}$. This explains the asymmetry of the signal and idler Schmidt modes in Fig. 4. The eigenvalues shown in Fig. 3 correspond to the geometric progression $q_{num}^{n}$, with $q_{num}=0.8903$, which is slightly higher than the value of $q$ given above by the Gaussian modeling. This and other slight differences between the numerical and the analytical solutions are caused by some of the four approximations we made.

The treatment of this section demonstrates the effectiveness of Gaussian modeling for the Schmidt modes, based on the complex Mehler’s formula.

V.3 Twin beams close to degeneracy

To illustrate the disappearance of the features peculiar to twin beams, now we consider a different angle between the crystal optical axis and the pump wave-vector of $\theta_{0}=$$$, for which the PDC is close to be frequency degenerate. The squeezing matrix is shown in Fig. 7. In this case the signal and the idler areas are not well separated and we expect deflections from the theory developed above.

The squeezing eigenvalues are shown in Fig. 8. Only the four first eigenvalues show multiplicity 2, predicted by Theorem 1.

The modal functions of the corresponding squeezing eigenmodes are shown in Fig. 9. Only the four first modal functions resemble concatenations of local Hermite-Gauss modes. At higher mode numbers the deflection of the transformation generator matrix from the form of Eq.(79) becomes significant and Theorem 1 loses its validity. Physically it means that the emission of two signal or two idler photons in one elementary act of photon-pair creation becomes significant as we approach the degenerate regime by varying the phase-matching conditions. We conjecture that the eigenvalues cease to be multiple at the mode number $k$ for which the Schmidt modal functions calculated separately for the signal and the idler beams start to overlap in the area around the origin. For the considered case it is $k=5$ and $k=6$, see the two bottom plots in Fig. 9, where the signal and the idler parts of the modal function start to deflect from the shape of the Hermite-Gauss function of the corresponding (second) order.

VI Conclusions

We have applied the formalism of Bloch-Messiah reduction to the parametric downconversion of light in the case where pulsed twin beams are generated. We have shown how the squeezing eigenvalues and the squeezing eigenmodes can be obtained in the case of moderate squeezing, where the solution of the wave equation is obtained in the first order of Magnus expansion. For this case we have proven a fundamental result: all the squeezing eigenvalues of twin beams have multiplicity at least two. As consequence, the modal functions of the squeezing eigenmodes are not unique and defined up to an orthogonal rotation in the space of two eigenmodes related to the same eigenvalue. We established two methods for avoiding the ambiguity in the definition of the squeezing eigenmodes: (i) reducing the Takagi factorization of the squeezing matrix to the spectral decomposition of an associated Hermitian matrix and (ii) tailoring the squeezing eigenmodes from the Schmidt modes of the signal and the idler beams obtained by SVD of the JSA matrix.

These general results have been illustrated by an example of twin beams discriminated by frequency. For the case of good separation of the signal and the idler spectra we have found the multiplicity two of all eigenvalues, predicted by the general theory. We have also shown that the Schmidt modal functions of the signal and the idler beams can be modeled with very good precision by replacing the JSA by a complex double-Gaussian function. The modal functions of the squeezing eigenmodes in this case are two chirped Hermite-Gauss functions in the signal and the idler frequency bands.

There are several important extensions of the present work. The considered example corresponds to rather narrowband twin beams discriminated by frequency. It would be interesting to apply the developed formalism to ultrabroadband fields with almost constant spectrum, for example, those generated in aperiodically poled crystals Horoshko and Kolobov (2013, 2017); Peřina (2018); Chekhova et al. (2018). In the narrowband regime a full treatment of 3D spatio-temporal modes is possible in the framework of complex analytic modeling developed here. It should reveal correlations between the spatial and temporal degrees of freedom of entangled twin beams. Another important direction of future research is related to twin beams discriminated by polarization or wave-vector direction. In the latter case two beams can be frequency-degenerate and their correlations can be converted to single-mode squeezing by means of interference on a beam splitter. This approach has high potential for producing highly multimode pulsed cluster states in a form of entangled frequency combs, similar to continuous-wave cluster states successfully produced in the last years Yokoyama et al. (2013).

Acknowledgments

D.B.H and M.I.K. thank G. Patera and T. Lipfert for helpful comments. F.A. is grateful to F. Corsi for useful discussions. This work was supported by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 665148 (QCUMbER). F.A. acknowledges financial support from the French National Research Agency Project No ANR-17-CE24-0035 VanQuTe.

Appendix A Proof of the complex Mehler’s formula

Here we prove a relation, which is equivalent to Eq. (V.2) and is obtained from it by letting $x^{\prime}=\tau_{1}x$, $y^{\prime}=\tau_{2}y$, $\eta^{\prime}=\eta/\sqrt{\mu\nu}$, $\xi^{\prime}=\xi/\sqrt{\mu\nu}$. Omitting the primes for simplicity we obtain the following identity

[TABLE]

which can be viewed as the Takagi factorization of a complex symmetric double-Gaussian kernel. Here $w=\sqrt{(1+\xi^{2})(1-\eta^{2})}$, $\zeta=\eta\xi/(2w)$, $p=\sqrt{1-q^{2}}$, $q=\sqrt{(1-q_{0})/(1+q_{0})}$, where $q_{0}=\sqrt{(1-\eta^{2})/(1+\xi^{2})}$. The angles $\theta$ and $\theta_{0}$ are found below. Note that the factor $\sqrt[4]{(1+\xi^{2})}$ in the right-hand side of Eq. (110) is the Hilbert-Schmidt norm of the kernel in its left-hand side.

Let us denote the kernel in the left-hand side of Eq. (110) by $K(x,y)$. This kernel corresponds to some integral operator $\mathcal{K}$. The standard method of finding the singular functions of this operator consists in solving the eigenvalue problems for the Hermitian operators $\mathcal{K}\mathcal{K}^{\dagger}$ and $\mathcal{K}^{\dagger}\mathcal{K}$. However, the phases of the singular functions cannot be determined in this way, since the eigenfunction is defined up to a unitary rotation (phase in one dimension), while the singular function in a Takagi factorization is defined up to an orthogonal rotation (sign in one dimension). For a full Takagi factorization we apply a more powerful method of generating function.

The generating function for the Hermite polynomials reads

[TABLE]

Now consider the following function

[TABLE]

where the last factor under the integral is the complement of the Hermite polynomial to the singular function in Eq. (110) (up to normalization). Differentiating the integrand $k$ times we obtain

[TABLE]

where $N_{k}=\sqrt{2^{k}k!\sqrt{\pi}}$ is the normalization constant of the Hermite-Gauss function.

On the other hand, Eq. (112) represents an integral of a complex Gaussian function of $y$, which can be taken analytically:

[TABLE]

where the complex numbers $q_{c}$ and $p_{c}$ are defined as

[TABLE]

Note that $p_{c}^{2}+q_{c}^{2}=1$. The last factor in the right-hand side of Eq. (114) is exactly $\Phi_{0}(x,q_{c}t)$. Differentiating it $k$ times we obtain

[TABLE]

Comparing Eq. (113) with Eq. (117) we conclude that the function $\Psi_{k}(y)=h_{k}(y)e^{i\zeta y^{2}+i(\theta_{0}+k\theta)/2}$, where $\theta_{0}=\arg(p_{c})$ and $\theta=\arg(q_{c})$, is the right-singular function of the kernel $K(x,y)$ with the singular value $|p_{c}q_{c}^{k}|$. It is not hard to find that $|q_{c}|=q$ and $|p_{c}|=p\sqrt[4]{(1+\xi^{2})}$. Since the considered kernel is symmetric, its left-singular function is $\Psi_{k}(x)$, which concludes the proof.