# Taking apart squeezed light

C. Drago<EMAIL_ADDRESS>J. E. Sipe
<EMAIL_ADDRESS>Department of Physics, University of Toronto, 60 St.
George Street, Toronto, Ontario, Canada, M5S 1A7

###### Abstract

We develop a formalism to describe squeezed light with large spectral-temporal
correlations. This description is valid in all regimes, but is especially
applicable in the long pulse to continuous-wave limit where the photon density
at any particular time is small, although the total number of photons can be
quite large. Our method relies on the Whittaker-Shannon interpolation formula
applied to the joint temporal amplitude of squeezed light, which allows us to
“take apart” the squeezed state. This provides a local description of the
state and its photon statistics, making the underlying physics more
transparent than does the use of the Schmidt decomposition. The formalism can
easily be extended to more exotic nonclassical states where a Schmidt
decomposition is not possible.

††preprint: APS/123-QED

##  I Introduction

Squeezed light is of interest for applications in quantum sensing and imaging
[1, 2], and as a resource for quantum computing [3]. For light propagating in
one direction in a quasi-1D structure, such as an optical fiber or a channel
waveguide in an integrated photonic structure [4], or even light propagating
in free space under the approximation that diffraction is negligible, a
squeezed state can be written as

$\displaystyle\left|\Psi\right\rangle=e^{\frac{\beta}{2}\int
d\omega_{1}d\omega_{2}\gamma(\omega_{1},\omega_{2})a^{\dagger}(\omega_{1})a^{\dagger}(\omega_{2})-\text{h.c.}}\left|\text{vac}\right\rangle,$
(I.1)

where for simplicity only one polarization and one transverse mode is
considered. We label the lowering operator at a frequency shifted by $\omega$
from a center reference frequency $\omega_{o}$ by $a(\omega)$ (see Appendix
A),

$\displaystyle[a(\omega_{1}),a^{\dagger}(\omega_{2})]=\delta(\omega_{1}-\omega_{2}),$
(I.2)

and $\left|\text{vac}\right\rangle$ is the vacuum state. Here
$\gamma(\omega_{1},\omega_{2})$ is the joint spectral amplitude,

$\displaystyle\int\left|\gamma(\omega_{1},\omega_{2})\right|^{2}d\omega_{1}d\omega_{2}=1,$
(I.3)

and $\beta$ is the squeezing amplitude; unless otherwise indicated, we take
integrals to range from $-\infty$ to $\infty$.

Often the properties of interest can be captured by simple functions of
frequencies, such as correlation functions of the form

$\displaystyle
G^{(1)}(\omega)=\left\langle\Psi|a^{\dagger}(\omega)a(\omega)|\Psi\right\rangle$
(I.4a) $\displaystyle
G^{(2)}(\omega_{1},\omega_{2})=\left\langle\Psi|a^{\dagger}(\omega_{1})a^{\dagger}(\omega_{2})a(\omega_{2})a(\omega_{1})|\Psi\right\rangle,$
(I.4b)

etc. For a pulse of light where $\left|\beta\right|\ll 1$, the state is only
slightly different from the vacuum state,

$\begin{split}\left|\Psi\right\rangle\approx|\text{vac}\rangle&+\frac{\beta}{2}\int
d\omega_{1}d\omega_{2}\gamma(\omega_{1},\omega_{2})a^{\dagger}(\omega_{1})a^{\dagger}(\omega_{2})\left|\text{vac}\right\rangle\end{split}$
(I.5)

where higher order terms in $\beta$ have been neglected, and there is only a
small probability amplitude for a two-photon state. Then

$\displaystyle
G^{(1)}(\omega)\rightarrow\left|\beta\right|^{2}\int\left|\gamma(\omega,\omega^{\prime})\right|^{2}d\omega^{\prime},$
(I.6a) $\displaystyle
G^{(2)}(\omega_{1},\omega_{2})\rightarrow\left|\beta\right|^{2}\left|\gamma(\omega_{1},\omega_{2})\right|^{2}.$
(I.6b)

In this paper we consider the evaluation of quantities such as these, even
when $\left|\beta\right|$ is not much less than one. A standard strategy in
such a situation is to decompose the joint spectral amplitude in terms of
Schmidt modes. If there is only one or a few Schmidt modes, as might occur for
squeezed light generated by a short pump pulse, the expressions for
correlation functions of the squeezed light in terms of Schmidt modes can be
easily evaluated even if $\left|\beta\right|$ is large, and they immediately
identify much of the physics. But large values of $\left|\beta\right|$ can
also arise for squeezed light generated by pump pulses that are very long, and
even if their intensities are very weak. Here, although the rate at which
pairs of photons are generated may be quite small, the total number of pairs
of photons generated diverges as the pump pulse approaches CW excitation, and
thus both $\left|\beta\right|$ and the Schmidt number diverge. Our goal is to
identify strategies that allow for the calculation of quantities such as
correlation functions to be done quickly for such states, and in a way that
makes the physics of the squeezed light clear.

We begin by introducing the temporal representation of the joint spectral
amplitude $\overline{\gamma}(t_{1},t_{2})$ and some of its general features in
Section II. Then in Section III we consider a natural first approach, which is
to use the Schmidt decomposition of $\overline{\gamma}(t_{1},t_{2})$ even if
the Schmidt number is very large. We find that this approach does not directly
make the physics of the state $\left|\Psi\right\rangle$ apparent, and this
motivates our search for other ways to “take apart” the joint amplitude of the
squeezed light that better elucidate the physics. In section IV we introduce a
new approach suggested by the time correlation functions of the light; it
works well if there is significant degeneracy in the amplitudes of the Schmidt
decomposition. In section V we generalize this, based on the Whittaker-Shannon
interpolation formula, in a scheme that is applicable even if there is no such
degeneracy. We argue that this new way of “taking apart” the joint amplitude
does make the physics of the state more apparent, and in section VI we compare
it to the Schmidt decomposition. Then using our formalism we provide a “local
decomposition” of the squeezed state, and demonstrate its use in calculations
in section VII. We give a discussion of the “strongly squeezed limit” within
our framework in section VIII, and end in section IX with a realistic example
of a joint spectral amplitude for a squeezed state generated in a ring
resonator structure. Our conclusions and suggestions for future work are
presented in section X.

##  II Joint temporal amplitude

Besides the expression (I.1) for a squeezed state that is based on an integral
over frequencies (or wavenumbers), it will be useful to have an expression
based on integrals over time (or position). For simplicity we assume group
velocity dispersion can be neglected and that light propagates with a velocity
$v$; then putting

$\displaystyle\overline{a}(t)=\int\frac{d\omega}{\sqrt{2\pi}}a(\omega)e^{-i\omega
t},$ (II.1)

and with

$\displaystyle\overline{\gamma}(t_{1},t_{2})\equiv\int\frac{d\omega_{1}d\omega_{2}}{2\pi}\gamma(\omega_{1},\omega_{2})e^{-i\omega_{1}t_{1}}e^{-i\omega_{2}t_{2}}$
(II.2)

identifying the “joint temporal amplitude,” we can write Eq. (I.1) as

$\displaystyle\left|\Psi\right\rangle=e^{\frac{\beta}{2}\int
dt_{1}dt_{2}\overline{\gamma}(t_{1},t_{2})\overline{a}^{\dagger}(t_{1})\overline{a}^{\dagger}(t_{2})-\text{h.c.}}\left|\text{vac}\right\rangle.$
(II.3)

We take $\gamma(\omega_{1},\omega_{2})$ and $\overline{\gamma}(t_{1},t_{2})$
to identify the spectral and temporal representations of a “joint amplitude”
and refer to their absolute squares $|\gamma(\omega_{1},\omega_{2})|^{2}$ and
$|\overline{\gamma}(t_{1},t_{2})|^{2}$ as the spectral and temporal
representations of a “joint intensity.” While due to time-ordering corrections
we would expect the joint amplitude to change as the pump intensity and thus
$\beta$ is increased [5], here we neglect such effects for simplicity, and
take the joint amplitude to be fixed when we consider varying $\beta$ below.

Of course, the ket $\left|\Psi\right\rangle$ (Eq. (I.1) or (II.3)) is a
Schrödinger ket at a particular time, say $t=0.$ The variables $t_{1},t_{2}$
can be thought of as surrogates for position, where $\overline{a}(t_{1})$ is
identified with the electric field at $z=-vt$1; equivalently, if the ket
$\left|\Psi\right\rangle$ were allowed to evolve in time, $\overline{a}(t)$
would identify the field operator at $z=0$ at time $t$ (see Appendix A).

Corresponding to the frequency correlation functions (Eq. (I.4)) we can also
introduce time dependent first- and second-order correlation functions [6]

$\displaystyle\overline{G}^{(1)}(t_{1},t_{2})=\left\langle\Psi|\overline{a}^{\dagger}(t_{1})\overline{a}(t_{2})|\Psi\right\rangle,$
(II.4a)
$\displaystyle\overline{G}^{(2)}(t_{1},t_{2})=\left\langle\Psi|\overline{a}^{\dagger}(t_{1})\overline{a}^{\dagger}(t_{2})\overline{a}(t_{2})\overline{a}(t_{1})|\Psi\right\rangle.$
(II.4b)

The “equal-time” first-order correlation function, given by
$\overline{G}^{(1)}(t)\equiv\overline{G}^{(1)}(t,t)$, is the “photon-density”
of the pulse of light, and is used to predict the counting rate of an ideal
photo-detector; indeed, if we integrate over all time, then

$N_{\text{pulse}}=\int\overline{G}^{(1)}(t)dt$ (II.5)

is the expected photon number in the pulse. The second-order correlation
function $\overline{G}^{(2)}(t_{1},t_{2})$ has a similar interpretation, and
is used to predict the probability of detection coincidences at the indicated
times.

Figure 1: Schematic of a joint intensity represented in time on the left and
frequency on the right. The horizontal width of the joint temporal (spectral)
amplitude is denoted by $\mathcal{T}_{p}$, ($\mathcal{B}_{c}$, the bandwidth)
is the effective pulse duration (bandwidth). The narrow horizontal width at
$t_{2}=0$ is denoted by $\mathcal{T}_{c}=1/\mathcal{B}_{c}$ and is the
coherence time of photon pairs.

We will primarily be interested in joint intensities of the form shown
schematically in Fig. 1, where there are large spectral-temporal correlations;
however, the formalism we introduce is valid for a general joint amplitude. We
characterize the joint intensity by two quantities, $\mathcal{T}_{p}$ and
$\mathcal{B}_{c}$, indicated in Fig. 1.

The quantity $\mathcal{T}_{p}$ is an effective measure of the duration of the
generated pulse of squeezed light; in practice this will depend on the pump
pulse duration if a nonresonant structure is used to generate the light, or by
the resonance time if a resonant structure is used.

The second quantity, $\mathcal{B}_{c}$, is an effective measure of the
bandwidth of generated photons; in practice this is set by phase-matching
constraints if a nonresonant structure is used to generate the light, or by
the resonator linewidth if a resonant structure is used. Associated with this
bandwidth we define a time $\mathcal{T}_{c}=1/\mathcal{B}_{c}$, which is on
the order of the narrow width of the joint temporal intensity; see Fig. 1. The
time $\mathcal{T}_{c}$ identifies the “coherence time,” the typical range of
$\left|t_{2}-t_{1}\right|$ over which
$\left|\overline{\gamma}(t_{1},t_{2})\right|^{2}$ is non-negligible.

##  III Schmidt modes

A natural approach to try to elucidate the physics of a squeezed state and
calculate the correlation functions given by Eq. (II.4) is to employ a Schmidt
decomposition of the joint amplitude, since the squeezed state can then be
written as a direct product of squeezed states associated with the supermodes
introduced via the Schmidt decomposition. That is the strategy we explore in
this section.

In writing Eq. I.1 for the squeezed state we took the origin of
$\gamma(\omega_{1},\omega_{2})$ to indicate the same frequency $\omega_{o}$
with respect to which both $\omega_{1}$ and $\omega_{2}$ are referenced – this
is the case of so-called “degenerate” squeezing – and here without loss of
generality the joint amplitude can be taken as symmetric in its variables
($\gamma(\omega_{2},\omega_{1})=\gamma(\omega_{1},\omega_{2})$, or
equivalently $\overline{\gamma}(t_{2},t_{1})=\overline{\gamma}(t_{1},t_{2})$).
From a Takagi factorization [7] we can then construct the Schmidt
decomposition,

$\begin{split}\gamma(\omega_{1},\omega_{2})&=\sum_{n}\sqrt{p_{n}}f_{n}(\omega_{1})f_{n}(\omega_{2}),\\\
\overline{\gamma}(t_{1},t_{2})&=\sum_{n}\sqrt{p_{n}}\overline{f}_{n}(t_{1})\overline{f}_{n}(t_{2}),\end{split}$
(III.1)

where the Schmidt weights $p_{n}\geq 0$ sum to unity; the elements of the sets
$\left\\{f_{n}(\omega)\right\\}$ and $\left\\{\overline{f}_{n}(t)\right\\}$ of
mutually orthogonal and normalized functions are related by

$\displaystyle\overline{f}_{n}(t)=\int\frac{d\omega}{\sqrt{2\pi}}f_{n}(\omega)e^{-i\omega
t}.$ (III.2)

As usual, we take the set $\left\\{f_{n}(\omega)\right\\}$ of elements with
$p_{n}\neq 0$ to be expanded to be a complete set of functions, assigning
$p_{n}=0$ to the added elements, and correspondingly for the
$\left\\{\overline{f}_{n}(t)\right\\}$.

The Schmidt number $K$ of the expansion (III.1), which characterizes the
effective number of spectral or temporal modes in the sum (III.1), is given by

$\displaystyle K=\left(\sum_{n}p_{n}^{2}\right)^{-1}.$ (III.3)

It serves as an effective measure of the correlation in the dependence of the
joint amplitude on its two variables; in the limit of a Schmidt number of
unity the joint amplitude is simply a product of the same function of each
variable, and there is no correlation at all in its dependence on those
variables.

We can introduce another measure of the correlation of the dependence of the
joint amplitude on its two variables by

$\mathcal{K}=\mathcal{T}_{p}\mathcal{B}_{c}=\frac{\mathcal{T}_{p}}{\mathcal{T}_{c}}.$
(III.4)

Clearly when $\mathcal{T}_{p}$ is much greater than $\mathcal{T}_{c}$ then the
dependence of the joint amplitude on its two variables is highly correlated
(see Fig. 1), and then $\mathcal{K}\gg 1$. Of course, at the moment we have
only defined $\mathcal{T}_{p}$ and $\mathcal{B}_{c}$ in a “rough-and-ready”
way, and such then is our definition of $\mathcal{K}$; we will make the
definition more precise below. Often similarly defined quantities are
introduced and referred to as the “time-bandwidth product.” So to avoid
confusion we henceforth refer to $\mathcal{K}$ as an “effective Schmidt
number.”

Confusions with the phrase “time-bandwidth product” can arise because of
different definitions used for “time” and “bandwidth.” For example, Fedorov et
al. calculate a time-bandwidth product by taking the width of the
“unconditional” (“single-particle”) and “conditional” (“coincidence”) widths
of the joint intensity; for a general double-Gaussian joint amplitude
(considered below), they show that the time-bandwidth product is exactly equal
to the Schmidt number [8, 9, 10]. Alternatively, Brecht and Silberhorn
calculate the time-bandwidth product taking the conditional and unconditional
widths of the “chronocyclic Wigner function,” which for their double-Gaussian
model is equal to the Schmidt number, even when a chirp is included [11]. The
agreement – either exact or approximate – between the “time-bandwidth product”
and the Schmidt number suggest a deeper connection between the two quantities.

However, there is an older and more rigorous meaning of the term “time-
bandwidth product” from classical information theory: For a one-dimensional
bandlimited signal, it is the number of orthogonal functions optimally
concentrated within a given timewidth needed to describe the signal [12], and
in this context is referred to as the “Shannon number” [13, 14].
Generalizations of the Shannon number exist for higher dimensional signals
where one calculates the time-bandwidth product using the bandwidth and
timewidth area, volume, etc. [13, 14, 15, 16]. Recent work has made the
connection between the Shannon number from classical information theory to the
Schmidt number in quantum information theory; see for example Pires et al.
[16] or Pors et al. [17] for both a theoretical and experimental
investigation.

In this spirit we define $\mathcal{K}$ more precisely by defining
$\mathcal{T}_{p}$ and $\mathcal{B}_{c}$ more precisely. We assume that those
quantities are chosen as small as possible but subject to the condition that,
to within the level of approximation adopted in calculations, they cover the
range of the joint amplitude in time and frequency respectively. This means
that the typical ways used to measure bandwidth or frequency – such as the
standard deviation, full-width-at-half-max etc. – are too narrow. In
particular, we choose $\mathcal{B}_{c}$ to be large enough that frequencies
larger than $2\pi\mathcal{B}_{c}/2$ can be completely neglected, and we can
treat the joint spectral amplitude as effectively bandwidth limited. For this
reason the effective Schmidt number $\mathcal{K}$ will generally be larger
than, but typically on the order of, other conventions used for the “time-
bandwidth product.”

In a later section we argue that generally the Schmidt number $K$ and
effective Schmidt number $\mathcal{K}$ satisfy the inequality

$K\leq\mathcal{K}.$ (III.5)

One might generally expect this to be true based on a physical argument:
Suppose we have squeezed light propagating with an arbitrary joint amplitude
with some Schmidt number $K$ and effective Schmidt number $\mathcal{K}$. Now
if the squeezed light is sent through a dispersive medium, the joint spectral
amplitude is multiplied by a complex but _separable_ phase that does _not_
change the Schmidt number. However, we know that in a dispersive medium the
bandwidth remains constant but the pulse duration broadens, and so
$\mathcal{K}$ will generally increase. Thus, in general one might indeed
expect that $K\leq\mathcal{K}$. Below we will discuss when the near exact
equality holds.

We will find that we can use the effective Schmidt number $\mathcal{K}$ to
introduce a “weak squeezing” regime characterized by the condition

$\frac{|\beta|}{\sqrt{\mathcal{K}}}\ll 1,$ (III.6)

and a “strong squeezing” regime characterized by the condition

$\frac{|\beta|}{\sqrt{\mathcal{K}}}\gg 1.$ (III.7)

If neither of these conditions are satisfied we refer to the squeezing as
“moderate.”

Now as a first example of the use of the Schmidt decomposition to evaluate the
correlation functions, consider a general normalized two-photon state,

$\displaystyle\left|\textup{II}\right\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\int
dt_{1}dt_{2}\overline{\gamma}(t_{1},t_{2})\overline{a}^{\dagger}(t_{1})\overline{a}^{\dagger}(t_{2})\left|\text{vac}\right\rangle$
(III.8)
$\displaystyle=\frac{1}{\sqrt{2}}\sum_{n}\sqrt{p_{n}}A_{n}^{\dagger}A_{n}^{\dagger}\left|\text{vac}\right\rangle,$

where we defined operators associated with the supermodes as

$\displaystyle A_{n}^{\dagger}\equiv\int
dt\overline{f}_{n}(t)\overline{a}^{\dagger}(t),$ (III.9)

and so

$\displaystyle\overline{a}(t)=\sum_{n}\overline{f}_{n}(t)A_{n}.$ (III.10)

To write the inverted form (Eq. (III.10)) we have taken the adjoint of Eq.
(III.9) and used the completeness relation of the set of functions
$\left\\{\overline{f}_{n}(t)\right\\}$. The set of operators
$\left\\{A_{n}\right\\}$ and their adjoints satisfy the usual harmonic
oscillator commutation relations. Just as the expressions (III.1) can be
written in time or frequency form, so the first of (III.8) and (III.9) can
also be written involving integrals over frequency of the corresponding
quantities. We find

$\begin{split}\left.\overline{G}^{(1)}(t)\right|_{\ket{\textup{II}}}&\equiv\left\langle\textup{II}|\overline{a}^{\dagger}(t)\overline{a}(t)|\textup{II}\right\rangle\\\
&=2\sum_{n}p_{n}\left|\overline{f}_{n}(t)\right|^{2},\end{split}$ (III.11)

and

$\begin{split}\left.\overline{G}^{(2)}(t_{1},t_{2})\right|_{\ket{\textup{II}}}&\equiv\left\langle\textup{II}|\overline{a}^{\dagger}(t_{1})\overline{a}^{\dagger}(t_{2})\overline{a}(t_{2})\overline{a}(t_{1})|\textup{II}\right\rangle\\\
&=2\left|\sum_{n}\sqrt{p_{n}}\overline{f}_{n}(t_{1})\overline{f}_{n}(t_{2})\right|^{2}.\end{split}$
(III.12)

In the first we have a contribution of $\left|\overline{f}_{n}(t)\right|^{2}$
from each Schmidt mode with a weight $p_{n}$, while in the second the
amplitudes associated with each of the Schmidt modes add; the “2” of course
arises because we have pairs of photons.

Moving to a squeezed state, in the limit $\left|\beta\right|\ll 1$ we can
write (II.3) as

$\displaystyle\left|\Psi\right\rangle\rightarrow\left|\text{vac}\right\rangle+\frac{\beta}{\sqrt{2}}\left|\textup{II}\right\rangle+.....,$
(III.13)

(cf. (I.5)), and we find

$\displaystyle\overline{G}^{(1)}(t)\rightarrow\sum_{n}\left|\beta_{n}\right|^{2}\left|\overline{f}_{n}(t)\right|^{2}=\left|\beta\right|^{2}\int\left|\overline{\gamma}(t,t^{\prime})\right|^{2}dt^{\prime}$
(III.14)
$\displaystyle\overline{G}^{(2)}(t_{1},t_{2})\rightarrow\left|\sum_{n}\beta_{n}\overline{f}_{n}(t_{1})\overline{f}_{n}(t_{2})\right|^{2}=\left|\beta\right|^{2}\left|\overline{\gamma}(t_{1},t_{2})\right|^{2},$

where $\beta_{n}=\beta\sqrt{p_{n}}$ (cf. (I.6)). Treating
$|\gamma(t_{1},t_{2})|^{2}$ as a normalized probability distribution, for
$N_{\text{pulse}}=|\beta|^{2}\ll 1$ we find $\overline{G}^{(1)}(t)$ is the
probability distribution reduced by integrating over the second time variable,
and $\overline{G}^{(2)}(t_{1},t_{2})$ is proportional to
$|\overline{\gamma}(t_{1},t_{2})|^{2}$ itself, the norm squared of the joint
temporal amplitude at the two corresponding times.

More generally, using the Schmidt decomposition (III.1) and the supermode
operators (III.9), which are all independent, we can write the squeezed ket
(II.3) as

$\displaystyle\left|\Psi\right\rangle=\bigotimes_{n}S_{n}\ket{\text{vac}}_{n},$
(III.15)

where $\ket{\text{vac}}_{n}$ is the vacuum state for the corresponding
supermode and

$\displaystyle
S_{n}=e^{\frac{\beta_{n}}{2}A_{n}^{\dagger}A_{n}^{\dagger}-\text{h.c.}}.$
(III.16)

With the standard result [18]

$\displaystyle
S_{n}^{\dagger}A_{n}S_{n}=c_{n}A_{n}+e^{i\theta}s_{n}A_{n}^{\dagger},$
(III.17)

where we have put $\beta=\left|\beta\right|e^{i\theta}$ and

$\displaystyle c_{n}\equiv\cosh\left(\left|\beta_{n}\right|\right),$ (III.18)
$\displaystyle s_{n}\equiv\sinh\left(\left|\beta_{n}\right|\right),$ (III.19)

using the expression (III.10) to write $\overline{a}(t)$ in terms of the
$\left\\{\overline{f}_{n}(t)\right\\}$ – and using (III.17) to evaluate
$\left\langle\Psi|A_{n}^{\dagger}A_{m}|\Psi\right\rangle$ and
$\left\langle\Psi|A_{n}^{\dagger}A_{m}^{\dagger}A_{p}A_{q}|\Psi\right\rangle$
– from (II.4) we have

$\displaystyle\overline{G}^{(1)}(t)=\sum_{n}\left|\overline{f}_{n}(t)\right|^{2}s_{n}^{2},$
(III.20)
$\displaystyle\overline{G}^{(2)}(t_{1},t_{2})=\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})+\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2}),$

where

$\displaystyle\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})=\left|\sum_{n}s_{n}c_{n}\overline{f}_{n}(t_{2})\overline{f}_{n}(t_{1})\right|^{2}$
(III.21a) $\displaystyle\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2})$
(III.21b) $\displaystyle\hskip
28.45274pt=\frac{1}{2}\sum_{n,m}\left|s_{n}s_{m}(\overline{f}_{n}(t_{1})\overline{f}_{m}(t_{2})+\overline{f}_{n}(t_{2})\overline{f}_{m}(t_{1}))\right|^{2},$

and the average photon number is

$\displaystyle N_{\text{pulse}}=\sum_{n}s_{n}^{2}.$ (III.22)

Of the two contributions to $\overline{G}^{(2)}(t_{1},t_{2})$, the “coherent”
term $\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})$ [19], which involves the
square of a sum, is the generalization to a squeezed state of the term
$\overline{G}^{(2)}(t_{1},t_{2})$ for the two-photon state (III.12), and the
only term that survives in $\overline{G}^{(2)}(t_{1},t_{2})$ in the limit
$\left|\beta\right|\ll 1$, cf. (III.14). The “incoherent” term
$\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2})$ [19] involves the sum of
squares, and will only be significant at larger values of
$\left|\beta\right|$. Note that the corresponding expressions for
$G^{(1)}(\omega)$ and $G^{(2)}(\omega_{1},\omega_{2})$ take the same form as
Eq. (III.20), with $\overline{f}_{n}$ replaced by $f_{n},$ $t_{1}$ by
$\omega_{1}$, etc.

###  III.1 Example 1: The double-Gaussian

Figure 2: For the double-Gaussian, from left to right we plot the: joint
temporal intensity divided by its maximum value with the axes normalized by
$\mathcal{T}_{p}^{\mathrm{DG}}$; joint spectral intensity divided by its
maximum value with the axes normalized by $2\pi\mathcal{B}_{c}^{\mathrm{DG}}$;
Schmidt amplitudes $p_{n}$ up to $n=34$. Figure 3: For the double-Gaussian,
from left to right we plot $\overline{G}^{(1)}(t)/\Phi$ and a few
contributions from different Schmidt modes in Eq. (III.20) with the horizontal
axis normalized by $\mathcal{T}_{p}^{\mathrm{DG}}$ for $\beta=0.1,5,$ and
$10$.

A simple model for the joint amplitude is a double-Gaussian function,

$\displaystyle\gamma(\omega_{1},\omega_{2})=\sqrt{\frac{1}{\pi\sigma_{p}\sigma_{c}}}e^{-\frac{(\omega_{1}-\omega_{2})^{2}}{4\sigma_{c}^{2}}}e^{-\frac{(\omega_{1}+\omega_{2})^{2}}{4\sigma_{p}^{2}}},$
(III.23a)
$\displaystyle\overline{\gamma}(t_{1},t_{2})=\sqrt{\frac{\sigma_{p}\sigma_{c}}{\pi}}e^{-\frac{\sigma_{c}^{2}(t_{1}-t_{2})^{2}}{4}}e^{-\frac{\sigma_{p}^{2}(t_{1}+t_{2})^{2}}{4}},$
(III.23b)

where $\sigma_{p}$ and $\sigma_{c}$ are the two parameters. The Schmidt modes
are harmonic oscillator wave functions. This can be seen by noting that the
reduced density operator of “particle 1” is equal to the density operator of a
harmonic oscillator in thermal equilibrium; its eigenfunctions are the Schmidt
modes, and they are obviously the harmonic oscillator wave functions. The
details can be worked out from this, or more mathematically from the Mehler
kernal [20, 21]. For $\sigma_{c}\geq\sigma_{p}$ , the Schmidt decompositions
are given by (III.1), with

$\displaystyle\gamma(\omega_{1},\omega_{2})=\sum_{n\geq
0}\sqrt{p_{n}}\mathcal{H}_{n}(\omega_{1})\mathcal{H}_{n}(\omega_{2}),$
(III.24a) $\displaystyle\overline{\gamma}(t_{1},t_{2})=\sum_{n\geq
0}\sqrt{p_{n}}\overline{\mathcal{H}}_{n}(t_{1})\overline{\mathcal{H}}_{n}(t_{2}),$
(III.24b)

and

$\displaystyle
p_{n}=\frac{4\sigma_{c}\sigma_{p}}{(\sigma_{c}+\sigma_{p})^{2}}\left(\frac{\sigma_{c}-\sigma_{p}}{\sigma_{c}+\sigma_{p}}\right)^{2n},$
(III.25)

with a Schmidt number

$\displaystyle K=\frac{\sigma_{c}^{2}+\sigma_{p}^{2}}{2\sigma_{c}\sigma_{p}}.$
(III.26)

Here $\overline{f}_{n}(t)=\overline{\mathcal{H}}_{n}(t)$, where

$\displaystyle\overline{\mathcal{H}}_{n}(t)=\frac{H_{n}(\frac{t}{t_{0}})e^{-t^{2}/(2t_{0}^{2})}}{\sqrt{2^{n}n!\pi^{1/2}t_{0}}},$
(III.27)

with $H_{n}(x)$ the Hermite polynomials, is the standard coordinate-
representation harmonic oscillator energy eigenfunction, but with $t$ playing
the role of $x$ and

$\displaystyle t_{0}\equiv\sqrt{\frac{1}{\sigma_{p}\sigma_{c}}}$ (III.28)

playing the role of a reference length $x_{0}$ that is often introduced [22];
$f_{n}(\omega)=\mathcal{H}_{n}(\omega),$ where

$\displaystyle\mathcal{H}_{n}(\omega)=(-i)^{n}\sqrt{\frac{t_{0}}{2^{n}n!\pi^{1/2}}}H_{n}(\omega
t_{0})e^{-\omega^{2}t_{0}^{2}/2}$ (III.29)

is the standard momentum-representation harmonic oscillator energy
eigenfunction [22], with $\omega$ playing the role of $p$.

In Fig. 2 we plot the joint intensities
$\left|\overline{\gamma}(t_{1},t_{2})\right|^{2}$ and
$\left|\gamma(\omega_{1},\omega_{2})\right|^{2}$ for
$\sigma_{c}/\sigma_{p}=50$, as well as the Schmidt weights $p_{n}$ as a
function of $n$. For $\sigma_{c}\gg\sigma_{p}$, the Schmidt number (III.26) is
approximately given by $K\approx\sigma_{c}/(2\sigma_{p})=25$, indeed we find
numerically that $K_{\mathrm{DG}}=25.01$. The joint intensities in Fig. 2 vary
over the widths $\mathcal{T}_{p}^{\mathrm{DG}}$ and
$\mathcal{B}_{c}^{\mathrm{DG}}$ which we set to be

$\mathcal{T}_{p}^{\mathrm{DG}}=\frac{a}{\sqrt{2}\sigma_{p}},\hskip
14.22636pt\mathcal{B}_{c}=\frac{a}{2\pi}\frac{\sigma_{c}}{\sqrt{2}},\hskip
14.22636pt\mathcal{T}_{c}^{\mathrm{DG}}=\frac{2\pi\sqrt{2}}{a\sigma_{c}},$
(III.30)

and choose $a=2\sqrt{2\pi}$ for convenience. This choice of $a$ is large
enough that the joint temporal and spectral amplitudes can essentially be
taken to be confined within the ranges of $\mathcal{T}_{p}^{\mathrm{DG}}$ and
$\mathcal{B}_{c}^{\mathrm{DG}}$ respectively, as can be gleaned from Fig. 2
and will in fact be confirmed by our calculations in later sections. Then the
effective Schmidt number is

$\mathcal{K}_{\mathrm{DG}}=\frac{a^{2}}{2\pi}\frac{\sigma_{c}}{2\sigma_{p}}\approx
4K_{\mathrm{DG}}=100.$ (III.31)

To illustrate the behaviour of the photon statistics for a large range of
photon numbers, we consider the three values of $\beta=0.1,5$, and $10$,
corresponding to $|\beta|/\sqrt{\mathcal{K}_{\mathrm{DG}}}=0.01,0.5,1$, and so
ranging from weak to moderate squeezing; we dedicate section VIII to the
discussion of the strongly squeezed limit.

In Fig. 3 we plot $\overline{G}^{(1)}(t)$ for the three values of $\beta$. The
expectation value of the number of photons is determined by (III.22), and for
$\beta=0.1,$ $5$, and $10$ we have respectively $N_{\text{pulse}}\approx
0.01,35$, and $383$; we take an effective photon flux (photons per unit time)
to be given by $\Phi=N_{\text{pulse}}/$$\mathcal{T}^{\mathrm{DG}}_{p}$. We
also plot the contribution to each $\overline{G}^{(1)}(t)$ from a number of
the Schmidt modes (see (III.20)). For each value of $\beta$, the photon
density at any particular time $t$ involves contributions from many Schmidt
modes, and we cannot associate it with one or even a few Schmidt modes.
Clearly as $\beta$ increases the contribution from the $n=0$ Schmidt mode
increases and the shape of the photon density narrows. This occurs because the
scaling of each contribution with $|\beta_{n}|$ is nonlinear and depends on
the quantities $c_{n}$ and $s_{n}$ (III.18); since for the double-Gaussian the
Schmidt amplitudes $p_{n}$ decrease as $n$ increases, $|\beta_{0}|$ has the
largest contribution. This behavior suggests that in the strongly squeezed
limit the photon statistics can be well approximated by the first few Schmidt
modes, a point to which we return in section VIII.

Figure 4: For the double-Gaussian, from left to right we plot
$\overline{G}^{(2)}(t_{1},t_{2})/\Phi^{2}$ (top) and the coherent and
incoherent contribution to $\overline{G}^{(2)}(t/2,-t/2)/\Phi^{2}$ (bottom)
with the axes normalized by $\mathcal{T}_{p}^{\mathrm{DG}}$ for $\beta=0.1,5,$
and $10$.

In Fig. 4 we turn to $\overline{G}^{(2)}(t_{1},t_{2})$ with again the three
values of $\beta$ considered above. In the top row we show this function for
the different values of $\beta$; at low $\beta$ the result is proportional to
the square of the absolute value of the joint temporal amplitude (see Fig. 2,
Eq. (III.14)), while for larger $\beta$ there are significant corrections to
this. In the bottom panel we plot $\overline{G}^{(2)}(t/2,-t/2)$, which
corresponds to moving along a diagonal that runs from the upper-left to the
lower-right of the plots in the first row; we give the coherent and incoherent
contributions separately.

The situation here is of course more complicated than that for
$\overline{G}^{(1)}(t)$, because the expression (III.20) for
$\overline{G}^{(2)}(t_{1},t_{2})$ is more complicated than a simple sum over
contributions from the individual Schmidt modes. But we see that at least in
the weak squeezing regime the coherent contribution to
$\overline{G}^{(2)}(t/2,-t/2)$ dominates, and it is nonzero only over a range
of $t$ much less than the range of $t$ over which the individual Schmidt modes
are nonzero; clearly the interference terms between the different Schmidt
modes in (III.20) play a critical role in the result for
$\overline{G}^{(2)}(t_{1},t_{2})$.

Further, the incoherent contribution at $\beta=5$ has a structure that
consists partly of a broad background and partly of a contribution that
mirrors the coherent contribution. Now note that the expression (III.21) for
$\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2})$ can also be very generally
written as

$\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2})=\overline{G}^{(1)}(t_{1})\overline{G}^{(1)}(t_{2})+|\overline{G}^{(1)}(t_{1},t_{2})|^{2},$
(III.32)

and the broad background in $\overline{G}_{\text{incoh}}^{(2)}(t/2,-t/2)$ can
be understood as arising from the first term on the right-hand-side. The
contribution that mirrors the coherent contribution can be understood as
arising from the second term, and in fact it is absent if we consider
nondegenerate squeezed light, where signal and idler frequencies are well
separated with different center frequencies [23]. The corresponding
contribution to $G^{(2)}_{\text{incoh}}(\omega_{1},\omega_{2})$ is often
discussed and referred to as the “autocorrelation” [24]. In any case, clearly
the different features of $\overline{G}^{(2)}(t/2,-t/2)$ certainly do not
follow in any simple way from the features of individual Schmidt modes, but
arise from the interference of many of these modes.

As $\beta$ increases, the range of $\overline{G}^{(2)}(t,t)$ narrows, as does
the range of the photon density. However, the coherent and incoherent
contribution along $\overline{G}^{(2)}(t/2,-t/2)$ broaden. Thus as the
squeezing parameter is increased, the photon statistics begin to look
uncorrelated. This behavior matches that of the photon density, in that as
$\beta$ increases fewer Schmidt modes are important in calculating the
correlation functions. We also note that the incoherent contribution is
approximately twice the coherent contribution, a point to which we return to
in section VIII.

While there are many features of interest here, we emphasize that both the
value of $\overline{G}^{(1)}(t)$ at a particular $t$, and that of
$\overline{G}^{(2)}(t_{1},t_{2})$ at a particular $t_{1}$ and $t_{2}$, receive
contributions from _many_ of the Schmidt modes. The same holds for the
corresponding functions $G^{(1)}(\omega)$ and
$G^{(2)}(\omega_{1},\omega_{2})$. The structure of the Schmidt modes
themselves, in and of itself, does not help us understand the structure of the
correlation functions.

Figure 5: For the sinc-hat, from left to right we plot the: joint temporal
intensity divided by its maximum value with the axes normalized by
$\mathcal{T}^{\mathrm{SH}}_{p}$; joint spectral intensity divided by its
maximum value with the axes normalized by $2\pi\mathcal{B}^{\mathrm{SH}}_{c}$;
Schmidt amplitudes $p_{n}$ up to $n=34$.  Figure 6: For the sinc-hat, from
left to right we plot $\overline{G}^{(1)}(t)/\Phi$ and a few contributions
from different Schmidt modes in Eq. (III.20) with the horizontal axis
normalized by $\mathcal{T}^{\mathrm{SH}}_{p}$ for $\beta=0.1,5,$ and $10$.

###  III.2 Example 2: The sinc-hat

This difference between the features of the Schmidt modes and the features of
$\overline{G}^{(1)}(t)$ and $\overline{G}^{(2)}(t_{1},t_{2})$ is not specific
to the double-Gaussian joint amplitude. Consider another form,

$\displaystyle\gamma(\omega_{1},\omega_{2})$
$\displaystyle=\alpha(\omega_{1}+\omega_{2})\phi\left(\frac{\omega_{1}-\omega_{2}}{2}\right),$
(III.33) $\displaystyle\overline{\gamma}(t_{1},t_{2})$
$\displaystyle=\overline{\alpha}\left(\frac{t_{1}+t_{2}}{2}\right)\overline{\phi}(t_{1}-t_{2}),$
(III.34)

where

$\displaystyle\alpha(\omega)=\frac{1}{\sqrt{\Omega_{p}}}\text{sinc}\left(\frac{\pi\omega}{\Omega_{p}}\right),$
(III.35) $\displaystyle\phi(\omega)=\frac{1}{\sqrt{\Omega_{c}}}\text{ for
}-\frac{\Omega_{c}}{2}\leq\omega\leq\frac{\Omega_{c}}{2},$ (III.36)
$\displaystyle=0,\text{ otherwise,}$
$\displaystyle\bar{\alpha}(t)=\frac{1}{\sqrt{T_{p}}}\text{ for
}-\frac{T_{p}}{2}\leq t\leq\frac{T_{p}}{2},$ (III.37) $\displaystyle=0,\text{
otherwise,}$
$\displaystyle\overline{\phi}(t)=\frac{1}{\sqrt{T_{c}}}\text{sinc}\left(\frac{\pi
t}{T_{c}}\right)$ (III.38)

with $T_{p}=2\pi/\Omega_{p}$ and $T_{c}=2\pi/\Omega_{c}$, and as usual
$\alpha(\omega)$ and $\overline{\alpha}(t),$ and $\phi(\omega)$ and
$\overline{\phi}(t)$, are Fourier transform pairs (cf. (III.2)). Since both
$\gamma(\omega_{1},\omega_{2})$ and $\overline{\gamma}(t_{1},t_{2})$ are
products of a “sinc” function and “top-hat” function we refer to this example
as the “sinc-hat” joint amplitude.

For the sinc-hat joint amplitude, the Schmidt modes must be found numerically.
In Fig. 5 we plot the joint intensities
$\left|\overline{\gamma}(t_{1},t_{2})\right|^{2}$ and
$\left|\gamma(\omega_{1},\omega_{2})\right|^{2}$, as well as the Schmidt
weights $p_{n}$ as a function of $n$, for $T_{p}/T_{c}=24$. We will see below
that for large $T_{p}/T_{c}$ we have $K\approx T_{p}/T_{c}$, and indeed here
we numerically find $K_{\text{SH}}=24.78$.

Evaluating $\overline{\gamma}(t_{1},t_{2})$ along the line $t_{1}=t_{2}$ those
variables range from $-T_{p}/2$ to $T_{p}/2$, and similarly along
$\omega_{1}=-\omega_{2}$, $\gamma(\omega_{1},\omega_{2})$ ranges from
$-\Omega_{c}/2$ to $\Omega_{c}/2$. Naïvely one would guess that we should set
$\mathcal{T}_{p}^{\text{SH}}\to T_{p}$ and
$2\pi\mathcal{B}_{c}^{\mathrm{SH}}\to\Omega_{c}$ or equivalently
$\mathcal{T}_{c}^{\text{SH}}\to T_{c}$, however, this is only the range of $t$
along the _diagonal_ (or anti-diagonal in frequency) and the joint amplitude
exists beyond it. In Appendix B we show that

$\mathcal{T}_{p}^{\text{SH}}=T_{p}+\frac{T_{c}}{2},$ (III.39)

and

$\mathcal{B}_{c}^{\mathrm{SH}}=\frac{\Omega_{c}}{2\pi}+\frac{1}{2}\frac{\Omega_{p}}{2\pi},\hskip
14.22636pt\mathcal{T}_{c}^{\mathrm{SH}}=\frac{T_{c}T_{p}}{T_{p}+T_{c}/2}.$
(III.40)

This leads to an effective Schmidt number

$\mathcal{K}_{\mathrm{SH}}=\frac{(T_{p}+\frac{T_{c}}{2})^{2}}{T_{p}T_{c}}=1+\frac{T_{p}}{T_{c}}+\frac{T_{c}}{4T_{p}},$
(III.41)

and for $T_{p}/T_{c}=24$, $\mathcal{K}_{\mathrm{SH}}\approx 25$ to very good
approximation. In section V we discuss this near equality.

In Fig. 6 we plot $\overline{G}^{(1)}(t)$ for the same three values of $\beta$
used in the example above, corresponding here to photon numbers
$N_{\text{pulse}}\approx 0.01,35$ and $335$; we take an effective photon flux
to be given by $\Phi=N_{\text{pulse}}/\mathcal{T}_{p}^{\mathrm{SH}}$. Then for
the three values of $\beta$ we have
$|\beta|/\sqrt{\mathcal{K}_{\mathrm{SH}}}\approx 0.02,1,2$, which again
corresponds to weak to moderate squeezing. In Fig. 6 we also include a few of
the contributions from the Schmidt modes; we see that typically those
contributions range over the whole duration of the pulse, analogous to what we
saw for the double-Gaussian example. In the top row of Fig. 7 we plot
$\overline{G}^{(2)}(t_{1},t_{2})$ for the indicated values of $|\beta|$, and
in the bottom row the coherent and incoherent contributions to
$\overline{G}^{(2)}(t/2,-t/2).$ Again, the range over which the individual
Schmidt modes extend is much larger than these contributions, and so they must
be understood as arising from a number of interfering Schmidt modes.

So just as for squeezed states described by the double-Gaussian joint
amplitude, the behavior of correlation functions of squeezed states described
by the sinc-hat joint amplitude cannot be linked in a simple way to the
behavior of the individual Schmidt modes. Quantitatively there are differences
between the correlation functions resulting from those two joint amplitudes:
The relative amplitudes of the Schmidt modes of a given $n$ in Fig. 6 (sinc-
hat joint amplitude) are roughly independent of $|\beta|$, while the relative
amplitudes in of those in Fig. 3 (double-Gaussian) are not, and the two parts
of the structure of $\overline{G}^{(2)}_{\text{incoh}}(t/2,-t/2)$ we noticed
for the double-Gaussian joint amplitude are even more pronounced, and persist
to larger $\beta$ than they did for that amplitude. These differences arise
because the Schmidt weights of the sinc-hat joint amplitude are nearly
identical in the $T_{p}/T_{c}\gg 1$ limit (see Fig. 5), and thus each Schmidt
mode contributes roughly equally to the resulting correlation functions even
in the large $\beta$ limit.

Figure 7: For the sinc-hat, from left to right we plot
$\overline{G}^{(2)}(t_{1},t_{2})/\Phi^{2}$ (top) and the coherent and
incoherent contribution to $\overline{G}^{(2)}(t/2,-t/2)/\Phi^{2}$ (bottom)
with the axes normalized by $\mathcal{T}^{\mathrm{SH}}_{p}$ for $\beta=0.1,5,$
and $10$.

To see how this near-degeneracy in the Schmidt weights arises, note that in
general the Schmidt modes $\overline{f}_{n}(t)$ of a joint temporal amplitude
$\overline{\gamma}(t_{1},t_{2})$ are eigenfunctions of the operator

$\displaystyle
M(t_{1},t_{2})\equiv\int\overline{\gamma}(t_{1},t)\overline{\gamma}^{*}(t,t_{2})dt,$
(III.42)

with eigenvalue $p_{n}$,

$\displaystyle\int
M(t_{1},t_{2})\overline{f}_{n}(t_{2})dt_{2}=p_{n}\overline{f}_{n}(t_{1}),$
(III.43)

which follows immediately from constructing $M(t_{1},t_{2})$ using (III.1).
Now for $T_{p}/T_{c}\gg 1$ we can approximate the joint temporal amplitude
(III.34) as

$\displaystyle\overline{\gamma}(t_{1},t_{2})$
$\displaystyle=\overline{\alpha}\left(\frac{t_{1}+t_{2}}{2}\right)\overline{\phi}(t_{1}-t_{2})$
(III.44)
$\displaystyle\approx\overline{\alpha}(t_{1})\overline{\phi}(t_{1}-t_{2}),$
(III.45)

and so

$\displaystyle M(t_{1},t_{2})$
$\displaystyle\approx\int\overline{\alpha}(t_{1})\overline{\alpha}(t_{2})\overline{\phi}(t_{1}-t)\overline{\phi}(t-t_{2})dt$
(III.46)
$\displaystyle=\sqrt{T_{c}}\overline{\alpha}(t_{1})\overline{\alpha}(t_{2})\overline{\phi}(t_{1}-t_{2})$
$\displaystyle=\frac{T_{c}}{T_{p}}\frac{\sin\left(\frac{\Omega_{c}}{2}(t_{1}-t_{2})\right)}{\pi(t_{1}-t_{2})},\text{
for }-\frac{T_{p}}{2}\leq t_{1},t_{2}\leq\frac{T_{p}}{2},$
$\displaystyle\text{ }=0,\text{ otherwise.}$

For specified $\Omega_{c}T_{p}/4$, the functions
$\overline{\psi}_{n}(t^{\prime})$ satisfying the eigenvalue equation

$\displaystyle\int\limits_{-T_{p}/2}^{T_{p}/2}\frac{\sin\left(\frac{\Omega_{c}}{2}(t-t^{\prime})\right)}{\pi(t-t^{\prime})}\overline{\psi}_{n}(t^{\prime})dt^{\prime}=\lambda_{n}\overline{\psi}_{n}(t),$
(III.47)

with the label $n=0,1,.....,$ are related to the angular prolate spheroidal
functions [12, 13, 14, 25, 26], and are defined with the normalization

$\displaystyle\int\limits_{-\infty}^{\infty}\left|\overline{\psi}_{n}(t)\right|^{2}dt=1.$
(III.48)

The $\lambda_{n}$ are close to unity for small $n$, and fall off quickly to
zero for $n>T_{p}/T_{c}\equiv K_{\text{app}};$ we approximate them as equal to
unity for $n<K_{\text{app}}$ and zero for $n\geq K_{\text{app}}$. So within
the approximation (III.46) we have

$\displaystyle\int
M(t_{1},t_{2})\overline{\psi}_{n}(t_{2})dt_{2}\rightarrow\frac{T_{c}}{T_{p}}\lambda_{n}\overline{\psi}_{n}(t_{1}),$
(III.49)

and comparing with (III.43) we can identify

$\displaystyle\overline{f}_{n}(t)\rightarrow\overline{\psi}_{n}(t),$ (III.50)
$\displaystyle
p_{n}\rightarrow\sqrt{\frac{T_{c}}{T_{p}}}=\frac{1}{\sqrt{K_{\text{app}}}},$

for $n<K_{\text{app}}$. That is, the Schmidt modes are approximately given by
the angular prolate spheroidal functions, and so

$\displaystyle\overline{\gamma}(t_{1},t_{2})\rightarrow\sum_{n=0}^{K_{\text{app}}-1}\sqrt{\frac{T_{c}}{T_{p}}}\overline{\psi}_{n}(t_{1})\overline{\psi}_{n}(t_{2}),$
(III.51)

exhibiting a huge degeneracy of Schmidt mode amplitudes, with a approximate
Schmidt number (III.3) $K_{\text{app}}=T_{p}/T_{c}\approx K_{\text{SH}}$, as
expected.

While (III.50,III.51) are only approximate (cf. Fig. 5), they do indicate that
in the limit $T_{p}/T_{c}\gg 1,$ which for the sinc-hat function we can
characterize as the “long pulse” limit, a large near-degeneracy of Schmidt
mode amplitudes can be expected. Like the exact Schmidt modes, the angular
prolate spheroidal functions range over the whole duration $T_{p}$ associated
with the joint temporal amplitude. However, were the degeneracy exact it would
imply that the Schmidt modes are not unique and that various superpositions of
them could be constructed. While this freedom is only approximate for near
degeneracy, we will see below that we can use it to “take apart” the joint
amplitude in a different way by constructing approximate Schmidt modes that
more explicitly reflect the properties of the light.

##  IV An approximate Schmidt decomposition

Focusing on the sinc-hat model (III.33,III.35) in the limit where $T_{p}\gg
T_{c}$ and $\mathcal{T}_{p}^{\mathrm{SH}}\to T_{p}$ and
$\mathcal{T}_{p}^{\mathrm{SH}}\to T_{c}$, note that while in the second row of
Fig. 7 we have plotted $\overline{G}^{(2)}(\overline{t}+t/2,\overline{t}-t/2)$
for $\overline{t}=0,$ we would expect such plots to be similar for values of
$\overline{t}\equiv(t_{1}+t_{2})/2$ ranging over the pulse duration,
especially in the limit of small $|\beta|$. In that limit
$\overline{G}^{(2)}(t_{1},t_{2})$ reflects the behavior of the joint temporal
amplitude itself (see (III.14)), and generally
$\overline{G}^{(2)}(\overline{t}+t/2,\overline{t}-t/2)$ will be nonzero for
$t$ ranging on a time scale of the order of $T_{c}$; in our example of Fig. 7
that is $T_{p}/24$. This suggests that if we want to capture the behavior of
the joint temporal amplitude as a function of $(t_{1}-t_{2})$ in each of the
terms of an approximate Schmidt decomposition, rather than just when they are
all used together, we should look for approximate Schmidt modes that vary over
a range of $T_{c}$. Since, roughly speaking, frequency components between
$-\Omega_{c}/2$ and $\Omega_{c}/2$ are then available, one such function can
easily be constructed by taking

$\begin{split}\overline{\eta}(t)&\equiv\frac{1}{\sqrt{\Omega_{c}}}\int\limits_{-\Omega_{c}/2}^{\Omega_{c}/2}\frac{d\omega}{\sqrt{2\pi}}e^{-i\omega
t}\\\ &=\frac{1}{\sqrt{T_{c}}}\text{sinc}\left(\frac{\pi
t}{T_{c}}\right),\end{split}$ (IV.1)

where the prefactor is chosen so the function is normalized (see (IV.3)
below). However, we need a set of such functions that are orthonormal to serve
as approximate Schmidt functions; the way to do that is to take the set of
functions

$\displaystyle\overline{\eta}_{n}(t)=\overline{\eta}(t-nT_{c}),$ (IV.2)

where $n$ is an integer, for then we have

$\displaystyle\int\overline{\eta}_{n}^{*}(t)\overline{\eta}_{m}(t)dt=\delta_{nm}.$
(IV.3)

Note that were time variables replaced by position variables, then the
$\left\\{\overline{\eta}_{n}(t)\right\\}$ would correspond to a set of Wannier
functions of the lowest band in a one-dimensional crystal of lattice spacing
corresponding to $T_{c}$, when the potential of the lattice is neglected
(“empty-lattice approximation”) [27]. We can then seek an approximate
expression $\overline{\gamma}_{\text{app}}(t_{1},t_{2})$ for the sinc-hat
joint temporal amplitude $\overline{\gamma}(t_{1},t_{2})$ of (III.33) by
writing

$\displaystyle\overline{\gamma}_{\text{app}}(t_{1},t_{2})=T_{c}\sum_{n}u(nT_{c})\overline{\eta}_{n}(t_{1})\overline{\eta}_{n}(t_{2}),$
(IV.4)

which clearly takes the form of a Schmidt decomposition, with $u(nT_{c})$
playing the role of an “envelope function”;
$\gamma_{\mathrm{app}}(t_{1},t_{2})$ is normalized as the exact function
(III.33) as long as

$\displaystyle T_{c}^{2}\sum_{n}\left|u(nT_{c})\right|^{2}=1.$ (IV.5)

The introduction of the functions $\overline{\eta}_{n}(t)$ allows us to work
with “pseudo-Schmidt” modes that are mutually orthogonal (like the real
Schmidt modes), but are localized in time and range over different center
times. We refer to the approximate Schmidt decomposition (IV.4) we construct
as the “pseudo-Schmidt decomposition.”

Figure 8: Plot of the approximate joint spectral intensity divided by its
maximum value with the axes normalized by $2\pi\mathcal{B}^{\mathrm{SH}}_{c}$.
The “false” contributions highlighted by the red dashed circles are due to the
periodicity of the function $\widehat{u}(\omega)$ with a period $T_{c}$; see
the discussion in the paragraph above Eq. (IV.14). Figure 9: From left to
right we plot the: sinc-hat joint temporal intensity divided by its maximum
value; approximate joint temporal intensity divided by its maximum value; two
contributions to the approximate joint temporal intensity divided by its
maximum value when $n=0$ and $n=5$ with the axes all normalized by
$\mathcal{T}^{\mathrm{SH}}_{p}$.

Taking the Fourier transform of the approximate joint temporal amplitude
(IV.4) gives the approximate joint spectral amplitude

$\displaystyle\gamma_{\text{app}}(\omega_{1},\omega_{2})$
$\displaystyle\equiv\int\frac{dt_{1}dt_{2}}{2\pi}\overline{\gamma}_{\text{app}}(t_{1},t_{2})e^{i\omega_{1}t_{1}}e^{i\omega_{2}t_{2}}$
(IV.6)
$\displaystyle=\frac{T_{c}^{2}}{2\pi}\widehat{u}(\omega_{1}+\omega_{2})s(\omega_{1})s\left(\omega_{2}\right),$

where

$\displaystyle\widehat{u}(\omega)\equiv\sum_{n}u(nT_{c})e^{in\omega T_{c}}$
(IV.7)

and

$\displaystyle s(\omega)$
$\displaystyle\equiv\frac{1}{\sqrt{T_{c}}}\int_{-\infty}^{\infty}\overline{\eta}(t)e^{i\omega
t}dt$ (IV.8) $\displaystyle=1,\text{ for
$-\frac{\Omega_{c}}{2}<\omega<\frac{\Omega_{c}}{2},$}$ $\displaystyle=0,\text{
otherwise.}$

We need to set $\widehat{u}(\omega),$ which satisfies
$\widehat{u}(\omega+mT_{c})=\widehat{u}(\omega)$ for any integer $m$. To
ensure that $\overline{\gamma}_{\text{app}}(t_{1},t_{2})$ is a good
approximation to $\overline{\gamma}(t_{1},t_{2})$ for the sinc-hat model, we
will want $u(nT_{c})$ to be independent of $n$ over an appropriate range.
Choosing $\mathcal{N}$ to be a large odd integer, we put

$\displaystyle u(nT_{c})$
$\displaystyle=\frac{1}{T_{c}\sqrt{\mathcal{N}}},\text{ for
$-\left(\frac{\mathcal{N}-1}{2}\right)\leq
n\leq\left(\frac{\mathcal{N}-1}{2}\right),$}$ $\displaystyle=0,\text{
otherwise.}$ (IV.9)

Then from (IV.7) this leads to

$\displaystyle\widehat{u}(\omega)=\frac{1}{T_{c}\sqrt{\mathcal{N}}}\frac{\sin\left(\frac{\mathcal{N}\omega
T_{c}}{2}\right)}{\sin\left(\frac{\omega T_{c}}{2}\right)}.$ (IV.10)

Now for the sinc-hat model we want the joint temporal amplitude
$\overline{\gamma}(t_{1},t_{2})$ to be nonvanishing for $t_{1},t_{2}$ varying
from $-T_{p}/2$ to $T_{p}/2$. From the approximate form
$\overline{\gamma}_{\text{app}}(t_{1},t_{2})$ that we are trying to construct
(IV.4), this implies that we should set

$\displaystyle T_{c}\left(\frac{\mathcal{N}-1}{2}\right)=\frac{T_{p}}{2},$
(IV.11)

or

$\displaystyle\mathcal{N}=1+\frac{T_{p}}{T_{c}}=\mathcal{K}_{\mathrm{SH}},$
(IV.12)

where the last equality holds to very good approximation when $T_{p}\gg
T_{c}$. For $T_{p}\gg T_{c}$, which is the limit we consider here, we can take
$\mathcal{N}$ to be either this or an odd integer close to it. The motivation
for calling $\mathcal{K}$ the “effective Schmidt number” is now clear, at
least for this joint amplitude; $\mathcal{K}_{\mathrm{SH}}$ is the effective
number of pseudo-Schmidt modes required for the decomposition. We can then
write (IV.10) as

$\displaystyle\widehat{u}(\omega)=\frac{1}{\sqrt{T_{c}T_{p}}}\frac{\sin\left(\frac{\omega}{2}(T_{c}+T_{c})\right)}{\sin\left(\frac{\omega}{2}T_{c}\right)}.$
(IV.13)

Note that with this choice we find
$\gamma_{\text{app}}(0,0)=\sqrt{T_{c}(T_{p}+T_{c})}/(2\pi)$, while for the
sinc-hat model we have exactly $\gamma(0,0)=\sqrt{T_{c}T_{p}}/(2\pi)$; so for
$T_{p}/T_{c}\gg 1$, $\gamma_{\text{app}}(0,0)$ is certainly equal to
$\gamma(0,0)$ to good approximation.

In Fig. 8 we plot
$\left|\gamma_{\text{app}}(\omega_{1},\omega_{2})\right|^{2}$ for
$T_{p}/T_{c}=24,$ taking $\mathcal{N}=\mathcal{K}_{\mathrm{SH}}=25$. Comparing
to the exact plot (middle plot in Fig. 5) we see that there is indeed
generally good agreement, with two main differences: First,
$\gamma_{\text{app}}(\omega_{1},\omega_{2})$ is only nonzero for
$(\omega_{1},\omega_{2})$ satisfying the bandwidth limiting conditions
$-\Omega_{c}/2\leq\omega_{1}\leq\Omega_{c}/2$ and
$-\Omega_{c}/2\leq\omega_{2}\leq\Omega_{c}/2$, while
$\gamma(\omega_{1},\omega_{2})$ extends beyond that; this is not so apparent
in Fig. 5, because outside the bandwidth limiting conditions the true
$\gamma(\omega_{1},\omega_{2})$ is very small. Second,
$\gamma_{\text{app}}(\omega_{1},\omega_{2})$ contains “false contributions”
near $(\omega_{1},\omega_{2})=(\Omega_{c}/2,\Omega_{c}/2)$ and
$(-\Omega_{c}/2,-\Omega_{c}/2$), see the contributions highlighted by the red
dashed circles in Fig. 8. This is related to the first difference, and arises
because the postulated approximate form (IV.4) of
$\overline{\gamma}_{\text{app}}(t_{1},t_{2})$ involves the function
$\widehat{u}(\omega)$ that is periodic in $\omega$ with a period $1/T_{c}$.

We now look at the correlation functions that follow for a squeezed state with
the approximate joint temporal amplitude identified here,

$\displaystyle\overline{\gamma}_{\text{app}}(t_{1},t_{2})={\displaystyle\sum_{n=-\frac{\mathcal{N}-1}{2}}^{\frac{\mathcal{N}-1}{2}}\sqrt{p_{n}}\overline{\eta}_{n}(t_{1})\overline{\eta}_{n}(t_{2}),}$
(IV.14)

where

$\displaystyle p_{n}$ $\displaystyle=\frac{1}{\mathcal{N}},\text{
for}\;-\left(\frac{\mathcal{N}-1}{2}\right)\leq
n\leq\left(\frac{\mathcal{N}-1}{2}\right)$
$\displaystyle=0\;\text{otherwise.}$ (IV.15)

From this point forward we will denote the sums over $n$ leaving the bounds
implicit. We plot $|\overline{\gamma}_{\text{app}}(t_{1},t_{2})|^{2}$ for the
sinc-hat model with $T_{p}/T_{c}=24$ in the middle diagram of Fig. 9,
repeating the exact $|\overline{\gamma}(t_{1},t_{2})|^{2}$ for this model in
the left-most diagram; we can see the general level of agreement that might be
expected from the results shown in Fig. 5 and Fig. 8. In the right-most plot
we show $\overline{\eta}_{n}(t_{1})\overline{\eta}_{n}(t_{2})$ for $n=0$ and
$n=5$; all such functions are of course well-localized, and from (IV.14) we
see that the contribution to $\overline{\gamma}_{\text{app}}(t_{1},t_{2})$
from each of these is their product with the (pseudo-) Schmidt weight $p_{n}$.

Using the general expressions (III.20,III.21) for $\overline{G}^{(1)}(t)$ and
$\overline{G}^{(2)}(t_{1},t_{2})$, we find that for our approximate model
(IV.14) we have

$\displaystyle\overline{G}^{(1)}(t)=s^{2}\sum_{n}\left|\overline{\eta}_{n}(t)\right|^{2},$
(IV.16)
$\displaystyle\overline{G}^{(2)}(t_{1},t_{2})=\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})+\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2}).$

where

$\displaystyle\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})$
$\displaystyle=s^{2}c^{2}\left|\sum_{p}\overline{\eta}_{p}(t_{2})\overline{\eta}_{p}(t_{1})\right|^{2},$
(IV.17) $\displaystyle\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2})$
$\displaystyle=\frac{1}{2}s^{4}\sum_{n,m}\left|\overline{\eta}_{n}(t_{1})\overline{\eta}_{m}(t_{2})+\overline{\eta}_{n}(t_{2})\overline{\eta}_{m}(t_{1})\right|^{2}.$

with

$s=\sinh\left(\frac{\left|\beta\right|}{\sqrt{\mathcal{N}}}\right),\hskip
14.22636ptc=\cosh\left(\frac{\left|\beta\right|}{\sqrt{\mathcal{N}}}\right).$
(IV.18)

These are of course the same formulas as for the exact Schmidt modes, except
that here $s$ and $c$ are the same for each pseudo-Schmidt mode. But since the
pseudo-Schmidt modes are localized over a time of order $T_{c}$, we see that
this decomposition of $\overline{G}^{(1)}(t)$ identifies contributions from
each pseudo-Schmidt mode that are localized in time windows much less than the
width of $\overline{G}^{(1)}(t)$, which is on the order of $T_{p}$. We show
some of these contributions in Fig. 10. Here the range over which $n$ varies
indicates the overall range of $\overline{G}^{(1)}(t)$ to very good
approximation, and at any given time $t$ the contributions to
$\overline{G}^{(1)}(t)$ come from at most a very few of the functions
$\overline{\eta}_{n}(t)$ with $nT_{c}$ close to $t$.

Figure 10: For the sinc-hat calculated from the pseudo-Schmidt decomposition,
from left to right we plot $\overline{G}^{(1)}(t)/\Phi$ and a few
contributions from different pseudo-Schmidt modes in Eq. (IV.16) with the
horizontal axis normalized by $\mathcal{T}^{\mathrm{SH}}_{p}$ for
$\beta=0.1,5,$ and $10$. Figure 11: For the sinc-hat calculated from the
pseudo-Schmidt decomposition, from left to right we plot
$\overline{G}^{(2)}(t_{1},t_{2})/\Phi^{2}$ (top) and the coherent and
incoherent contribution to $\overline{G}^{(2)}(t/2,-t/2)/\Phi^{2}$ (bottom)
with the axes normalized by $\mathcal{T}^{\mathrm{SH}}_{p}$ for $\beta=0.1,5,$
and $10$.

At least within the approximate pseudo-Schmidt decomposition we can now
justify the use of the condition (III.6) to identify the weak squeezing
regime. Since in our example here $\mathcal{N}=\mathcal{K}_{\text{SH}}$, the
weak squeezing limit can be written as

$\frac{|\beta|}{\sqrt{\mathcal{N}}}\ll 1,$ (IV.19)

the expectation value of the number of photons in any pseudo-Schmidt mode is
given by

$\displaystyle
N_{\text{mode}}=s^{2}=\sinh^{2}\left(\frac{\left|\beta\right|}{\sqrt{\mathcal{N}}}\right),$
(IV.20)

while the expectation value of the number of photons in the full pulse is

$\displaystyle
N_{\text{pulse}}=\mathcal{N}\sinh^{2}\left(\frac{\left|\beta\right|}{\sqrt{\mathcal{N}}}\right).$
(IV.21)

So in the weak squeezing limit (IV.19) we have $N_{\text{mode}}\ll 1$; the
expected number of photons in any pseudo-Schmidt mode is much less than unity.
However, the expected number of photons in the pulse, $N_{\text{pulse}}$ can
be arbitrarily large; indeed, in the limit of weak squeezing we have
$N_{\text{pulse}}\rightarrow|\beta|^{2}.$ Very generally, regardless of the
level of squeezing, we can understand the CW limit as
$|\beta|^{2}\rightarrow\infty$ and $\mathcal{N}\rightarrow\infty$ such that
$|\beta|^{2}/\mathcal{N}$ is a constant.

In Fig. 11 we plot the approximate $\overline{G}^{(2)}(t_{1},t_{2})$ given by
(IV.16,IV.17), and compare with the exact $\overline{G}^{(2)}(t_{1},t_{2})$ in
Fig. 7 for the sinc-hat model; we see very good agreement. Importantly, the
decomposition (IV.16,IV.17) of $\overline{G}^{(2)}(t_{1},t_{2})$ into pseudo-
Schmidt modes immediately illustrates the behavior of that function in a way
that the decomposition into the exact Schmidt modes does not.

Considering the first plot when $\beta=0.1$ ($|\beta|/\sqrt{\mathcal{N}}$ =
0.02), the dominance of $\overline{G}^{(2)}_{\text{coh}}(t/2,-t/2)$ over
$\overline{G}^{(2)}_{\text{incoh}}(t/2,-t/2)$ follows immediately from the
prefactors $s^{2}c^{2}$ and $s^{4}$ in their expressions (IV.17). And in the
sum over $p$ in the expression for
$\overline{G}^{(2)}_{\text{coh}}(t/2,-t/2)$, of the terms
$\overline{\eta}_{p}(t/2)\overline{\eta}_{p}(-t/2)$ both
$\overline{\eta}_{p}(t/2)$ and $\overline{\eta}_{p}(-t/2)$ must be significant
for a contribution to be made, which requires $p$ to be close to $0$ and
$|t|\lesssim T_{c}$. The same behavior extends for the larger values of
$\beta$ as in the exact sinc-hat Schmidt decomposition. Thus the largest
contribution to $\overline{G}^{(2)}(t/2,-t/2)$ when $|t|\lesssim T_{c}$ comes
from only a few terms in the pseudo-Schmidt decomposition, while it involves
many terms in the Schmidt decomposition, and their interference.

The expression for $\overline{G}^{(2)}_{\text{incoh}}(t/2,-t/2)$, which in the
pseudo-Schmidt decomposition (IV.17) involves sums over two indices $n$ and
$m$, contains two types of contributions. In the first, with $m=n$, we get
terms that will only be significant if $m=n$ is close to $0$ and $|t|\lesssim
T_{c}$, as in the expression for $\overline{G}^{(2)}_{\text{coh}}(t/2,-t/2)$;
this gives the contribution to $\overline{G}^{(2)}_{\text{incoh}}(t/2,-t/2)$
that mirrors the form of $\overline{G}^{(2)}_{\text{coh}}(t/2,-t/2)$. This
contribution to $\overline{G}^{(2)}_{\text{incoh}}(t/2,-t/2)$, and the term
$\overline{G}^{(2)}_{\text{coh}}(t/2,-t/2)$, can thus be seen to arise from
pairs of photons, each photon in a pair associated with the same pseudo-
Schmidt mode. But the terms with $m\neq n$ can give contributions for $t$ on
the order of $T_{p}$; they give rise to the broad background, which can be
understood as arising from pairs of photons, with the photons in a pair
associated with different pseudo-Schmidt modes.

In a similar way one can understand the behavior of
$\overline{G}^{(2)}_{\text{coh}}(\overline{t}+t/2,\overline{t}-t/2)$ and
$\overline{G}^{(2)}_{\text{incoh}}(\overline{t}+t/2,\overline{t}-t/2)$ for
$\overline{t}\neq 0$. Unlike in the decomposition of the joint temporal
amplitude in terms of Schmidt modes, the decomposition in terms of pseudo-
Schmidt modes immediately reveals the structure of those correlation
functions.

And in fact, we can derive an analytic expression for the approximate
(IV.16,IV.17) $\overline{G}^{(1)}(t)$ and $\overline{G}^{(2)}(t_{1},t_{2})$ in
the CW limit. Noting that

$\displaystyle\sum_{n=-\infty}^{\infty}\overline{\eta}_{n}(t_{2})\overline{\eta}_{n}(t_{1})=\frac{1}{T_{c}}\text{sinc}\left(\frac{\pi\>\Delta
t}{T_{c}}\right),$ (IV.22)

where $\Delta t\equiv t_{2}-t_{1}$, for $t_{1}$ and $t_{2}$ in the center of a
pulse of duration $T_{p}$, when $T_{p}\rightarrow\infty$ we can write

$\displaystyle\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})\rightarrow\frac{s^{2}c^{2}}{T_{c}^{2}}\text{sinc}^{2}\left(\frac{\pi\>\Delta
t}{T_{c}}\right),$ (IV.23)

while

$\displaystyle\overline{G}^{(1)}(t)\rightarrow\frac{s^{2}}{T_{c}}.$ (IV.24)

Finally, noting that

$\begin{split}\overline{G}^{(1)}(t_{1},t_{2})&=s^{2}\sum_{n}\overline{\eta}_{n}(t_{2})\overline{\eta}_{n}(t_{1})\\\
&\rightarrow\frac{s^{2}}{T_{c}}\text{sinc}\left(\frac{\pi\>\Delta
t}{T_{c}}\right),\end{split}$ (IV.25)

then with the alternate form of $\overline{G}_{\text{incoh}}^{(2)}$ (see Eq.
(III.32)) we have

$\displaystyle\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2})\rightarrow\frac{s^{4}}{T_{c}^{2}}\left(1+\text{sinc}^{2}\left(\frac{\pi\>\Delta
t}{T_{c}}\right)\right),$ (IV.26)

and so

$\displaystyle\overline{G}^{(2)}(t_{1},t_{2})$
$\displaystyle\rightarrow\overline{G}^{(2)}(\Delta
t)=\frac{s^{2}c^{2}}{T_{c}^{2}}\text{sinc}^{2}\left(\frac{\pi\>\Delta
t}{T_{c}}\right)$ (IV.27) $\displaystyle\hskip
42.67912pt+\frac{s^{4}}{T_{c}^{2}}\left(1+\text{sinc}^{2}\left(\frac{\pi\>\Delta
t}{T_{c}}\right)\right).$

We plot this in Fig. 12 (cf. Fig. 7 and 11.)

From the above discussion, we are motivated to think of
$\overline{G}^{(2)}(\overline{t}+t/2,\overline{t}-t/2)$ for times $|t|\lesssim
T_{c}$ when $n=m$ as arising from one pseudo-Schmidt mode at a time, as it
were, on a “mode-by-mode” basis. Since the pseudo-Schmidt decomposition is
valid only for joint amplitudes such as the sinc-hat joint amplitude, where
there is a high degeneracy in the Schmidt weights, we give this calculation in
Appendix F.

Figure 12: Using the analytical result (IV.27), from left to right we plot
$\overline{G}^{(2)}(t_{1},t_{2})/\Phi^{2}$ (top) and the coherent and
incoherent contribution to $\overline{G}^{(2)}(t/2,-t/2)/\Phi^{2}$ (bottom)
with the axes normalized by $\mathcal{T}^{\mathrm{SH}}_{p}$ for $\beta=0.1,5,$
and $10$.

##  V The Whittaker-Shannon decomposition

The ability to construct approximate pseudo-Schmidt modes for the sinc-hat
joint amplitude that are “localized” in time ($\overline{\eta}(t)$, see
(IV.1)) relied on the near-degeneracy of the Schmidt modes. But this cannot be
generally expected; see, e.g., Fig. 2 for the Schmidt amplitudes of the
double-Gaussian joint amplitude, where there is no such near-degeneracy.
Nonetheless, the recognition used above that only a finite frequency range –
there from $-\Omega_{c}/2$ to $\Omega_{c}/2$ – is important can be employed to
construct an extension of a Schmidt decomposition more generally; it shares
the feature of the pseudo-Schmidt decomposition above that the functions
involved are localized in time.

We consider a joint spectral amplitude $\gamma(\omega_{1},\omega_{2})$ to be
bandwidth limited, in that at least approximately it can be taken to be
nonzero only for $-\Omega/2\leq\omega_{1},\omega_{2}\leq\Omega/2$, where
$\Omega$ is a positive frequency. Then we can construct a Whittaker-Shannon
decomposition of the joint amplitude based on the Whittaker-Shannon Sampling
Theorem [28, 29, 30]. That theorem states that for a function
$\overline{g}(t)$ that is bandwidth limited in the sense defined above – i.e.,
involving only constituent frequencies in the range
$-\Omega/2\leq\omega\leq\Omega/2$ – we can write

$\displaystyle\overline{g}(t)=\sum_{n}\overline{g}(n\tau)\text{sinc}\left(\frac{(t-n\tau)\pi}{\tau}\right),$
(V.1)

where $n$ ranges over the integers and

$\displaystyle\tau=\frac{2\pi}{\Omega}.$ (V.2)

Defining

$\displaystyle\overline{\chi}(t)=\frac{1}{\sqrt{\tau}}\text{sinc}\left(\frac{\pi
t}{\tau}\right)$ (V.3)

and putting

$\displaystyle\overline{\chi}_{n}(t)=\overline{\chi}(t-n\tau),$ (V.4)

we have

$\displaystyle\int\overline{\chi}_{n}^{*}(t)\overline{\chi}_{m}(t)dt=\delta_{nm},$
(V.5)

and we can write (V.1) as

$\displaystyle\overline{g}(t)=\sqrt{\tau}\sum_{n}\overline{g}(n\tau)\overline{\chi}_{n}(t).$
(V.6)

The Fourier transform of $\overline{g}(t)$ is

$\displaystyle g(\omega)=\int\frac{dt}{\sqrt{2\pi}}\overline{g}(t)e^{i\omega
t}=\sqrt{\tau}\sum_{n}\overline{g}(n\tau)\chi_{n}(\omega),$ (V.7)

where

$\displaystyle\chi_{n}(\omega)$ $\displaystyle=\frac{e^{i\omega
n\tau}}{\sqrt{\Omega}}\text{ for
$-\frac{\Omega}{2}\leq\omega\leq\frac{\Omega}{2}$}$ (V.8)
$\displaystyle=0\text{ otherwise,}$ (V.9)

is the Fourier transform of $\overline{\chi}_{n}(t)$; we have

$\displaystyle\int\chi_{n}^{*}(\omega)\chi_{m}(\omega)d\omega=\delta_{nm}.$
(V.10)

For frequencies within the band limited region the set of functions
$\\{\chi_{n}(\omega)\\}$ and $\\{\chi_{n}(t)\\}$ form an approximately
complete set of functions and we refer to them as the “Whittaker-Shannon
modes.”

The “Whittaker-Shannon decomposition” of the joint temporal amplitude
$\overline{\gamma}(t_{1},t_{2})$, and its Fourier transform
$\gamma(\omega_{1},\omega_{2})$, follow immediately from this if
$\gamma(\omega_{1},\omega_{2})$ is bandwidth limited, with frequency bandwidth
$\Omega.$ Using the sampling theorem for both variables, we have

$\displaystyle\overline{\gamma}(t_{1},t_{2})$
$\displaystyle=\tau\sum_{n,m}\overline{\gamma}(n\tau,m\tau)\overline{\chi}_{n}(t_{1})\overline{\chi}_{m}(t_{2}),$
(V.11) $\displaystyle\gamma(\omega_{1},\omega_{2})$
$\displaystyle=\tau\sum_{n,m}\overline{\gamma}(n\tau,m\tau)\chi_{n}(\omega_{1})\chi_{m}(\omega_{2}),$

Unlike a Schmidt (or pseudo-Schmidt) decomposition, this involves a double
sum. But the Whittaker-Shannon supermodes involved, with raising operators
given by

$\displaystyle B_{n}^{\dagger}=\int
d\omega\chi_{n}(\omega)a^{\dagger}(\omega)=\int
dt\overline{\chi}_{n}(t)\overline{a}^{\dagger}(t),$ (V.12)

are associated with functions that are localized in time
($\overline{\chi}_{n}(t))$ and can be inverted so that

$\overline{a}(t)=\sum_{n}\overline{\chi}_{n}(t)B_{n}.$ (V.13)

It will be useful to define

$\beta_{nm}=\beta\tau\overline{\gamma}(n\tau,m\tau),$ (V.14)

which is generally complex but symmetric, $\beta_{nm}=\beta_{mn}$; then we can
write the squeezed state $\left|\Psi\right\rangle$ (I.1,II.3) as

$\displaystyle\left|\Psi\right\rangle=e^{\frac{1}{2}\sum\limits_{n,m}\beta_{nm}B_{n}^{\dagger}B_{m}^{\dagger}-\text{h.c.}}\ket{\text{vac}}.$
(V.15)

Although non-diagonal terms of $\beta_{nm}$ will be important, since
$\overline{\gamma}(t_{1},t_{2})$ is only significant for
$\left|t_{1}-t_{2}\right|$ on the order of the coherence time, we can expect
$\beta_{nm}$ to be significant only for $m$ “reasonably close” to $n$; we will
examine this in more detail below. Note that our examples of the double-
Gaussian and sinc-hat joint spectral amplitudes illustrate that imposing the
approximation that those amplitudes are bandwidth limited requires a choice of
$\Omega\geq 2\pi\mathcal{B}_{c}$ that corresponds to
$\tau\leq\mathcal{T}_{c}$, that is, $\tau$ is smaller than, but on the order
of, the coherence time.

As a first example, we apply the Whittaker-Shannon decomposition to the sinc-
hat joint amplitude (III.33,III.35), and show how the approximate pseudo-
Schmidt decomposition (IV.14,IV) arises as the limiting case when $T_{p}\gg
T_{c}$. To apply the Whittaker-Shannon decomposition we must choose a
bandlimit, and for the sinc-hat example the joint spectral amplitude ranges
mostly within the “box” set by the width $\Omega_{c}$; a natural choice for
the bandlimit $\Omega$ is then to set $\Omega\to\Omega_{c}$ and it follows
that $\tau\to T_{c}$. However, the sinc-hat joint spectral amplitude will
never fit exactly inside a box of width $\Omega_{c}$ because it is only
exactly bandlimited in the diagonal direction set by the range of
$\phi((\omega_{1}-\omega_{2})/2)$, see Eq. (III.33), (III.35), so the corners
will always exist outside the box, see the discussion surrounding Eq. (III.40)
and Appendix B. When we work in the limit that $T_{p}\gg T_{c}$, where the
contributions that exist outside the boundary of the box are very small, then
to good approximation we can treat the sinc-hat joint spectral amplitude as
bandlimited with bandwidth $\Omega=\Omega_{c}$ ($\tau=T_{c}$).

With this choice of $\Omega$ and the sinc-hat joint temporal amplitude (Eq.
(III.33)), we evaluate

$\begin{split}\tau\>\overline{\gamma}(n\tau,m\tau)&=\delta_{nm}\sqrt{\frac{T_{c}}{T_{p}}},\text{for}-\frac{T_{p}}{2}\leq
nT_{c}\leq\frac{T_{p}}{2}\\\ &=0,\text{otherwise},\end{split}$ (V.16)

where we have used $\text{sinc}((n-m)\pi)=\delta_{nm}$. We find that the range
of Whittaker-Shannon modes is then diagonal, and only ranges over the pulse
duration set by $\overline{\alpha}(nT_{c})$; identifying that
$T_{p}/T_{c}=\mathcal{N}-1$ (IV.12) and $p_{n}$ (IV) as in section IV, it
immediately follows that

$\displaystyle\overline{\gamma}(t_{1},t_{2})\rightarrow\sum_{n=-\left(\frac{\mathcal{N}-1}{2}\right)}^{\left(\frac{\mathcal{N}-1}{2}\right)}\sqrt{p_{n}}\overline{\chi}_{n}(t_{1})\overline{\chi}_{n}(t_{2}),$
(V.17)

where $\overline{\chi}_{n}(t)=\overline{\eta}_{n}(t)$.

For the sinc-hat model (see Fig. 5), the reduction of the double sum to its
diagonal contributions is only possible because the joint spectrum is exactly
bandlimited in the $(\omega_{1}-\omega_{2})$ direction (see (III.33) and
(III.35)); however, the resulting Whittaker-Shannon decomposition will only be
accurate as a single sum in the $T_{p}\gg T_{c}$ limit for which the
contributions that exist outside the box become very small and can safely be
neglected.

In general, this reduction will not be possible and we will need to include
some – but typically only a few – off diagonal contributions to properly
approximate the joint amplitude with the Whittaker-Shannon decomposition. As
we will see below, the single sum and product state that follows if only
diagonal terms are included (see Eq. (III.15)) is much easier and more
intuitive to work with, so we always want to be as close to that limit as we
can be. This means that for a not-exactly-bandlimited joint spectral
amplitude, the choice of a bandwidth $\Omega=2\pi/\tau$ always involves a
trade-off: We want $\tau$ as large as possible so that each
$\overline{\gamma}(n\tau,n\tau)$ covers most of the $|t_{1}-t_{2}|$ behavior;
however, if we choose $\tau$ too large then $\Omega$ becomes too small to even
approximately cover the bandwidth of the joint spectral amplitude.

To see how this plays out in practice, consider the double-Gaussian joint
amplitude (III.23). Although it is not strictly bandlimited, if we set
$\Omega/2\pi=\mathcal{B}_{c}^{\mathrm{DG}}=a\sigma_{c}/(2\pi\sqrt{2})$ (see
Eq. III.30), to very good approximation we can neglect the high frequency
contributions for a reasonable choice of $a$. Then
$\tau=2\pi\sqrt{2}/(a\sigma_{c})$, and

$\overline{\gamma}(n\tau,m\tau)\propto
e^{-\frac{\sigma_{c}^{2}\tau^{2}(n-m)^{2}}{4}}$ (V.18)

sets the off diagonal range; for our choice of $\tau$ we have
$\sigma_{c}^{2}\tau^{2}=8\pi^{2}/a^{2}$ and as $a$ increases the bandlimit
gets larger, but so does the range over which $|n-m|$ is significant. Setting
$a=2\sqrt{2\pi}$ as in Section III, we have $\sigma_{c}^{2}\tau^{2}/4=\pi$ and
the $1/e$ drop-off in Eq. (V.18) occurs when $|n-m|=1/\sqrt{\pi}\approx 0.5$.
In this example, and more generally for only approximately band limited joint
spectral amplitudes, one could investigate optimizing $\tau$ constrained by a
specified error tolerance on the Whittaker-Shannon interpolation.

Clearly the Whittaker-Shannon decomposition depends on the two index parameter
$\beta_{nm}$, but to make a comparison to the pseudo-Schmidt decomposition we
focus on $\mathring{\beta}$, which we define to be the value of $\beta_{nm}$
at the $n$ and $m$ for which $|\overline{\gamma}(n\tau,m\tau)|$ takes its
maximum value. Denoting the value of $|\overline{\gamma}(n\tau,m\tau)|$ at
this $n$ and $m$ by $\overline{\gamma}_{\text{max}}$, we have

$\mathring{\beta}=\tau\beta\overline{\gamma}_{\text{max}},$ (V.19)

and

$\beta_{nm}=\mathring{\beta}r_{nm},$ (V.20)

with $r_{nm}=\overline{\gamma}(n\tau,m\tau)/\overline{\gamma}_{\text{max}}$.
The range of $|n-m|$ over which $r_{nm}$ is significant identifies the range
over which $\beta_{nm}$ varies.

Again using the double-Gaussian joint amplitude as an example, which achieves
its maximum at $t_{1}=t_{2}=0$, we find

$|\mathring{\beta}|=\sqrt{2}|\beta|\frac{\tau}{\sqrt{\mathcal{T}^{\text{DG}}_{p}\mathcal{T}^{\text{DG}}_{c}}},$
(V.21)

where we have used
$\mathcal{T}^{\text{DG}}_{p}\mathcal{T}^{\text{DG}}_{c}=2\pi/\sigma_{c}\sigma_{p}$.
From the discussion above, we set $\tau=\mathcal{T}_{c}^{\text{DG}}$, and
obtain

$|\mathring{\beta}|=\sqrt{2}|\beta|\sqrt{\frac{\mathcal{T}_{c}^{\text{DG}}}{\mathcal{T}_{p}^{\text{DG}}}}=\sqrt{2}\frac{|\beta|}{\sqrt{\mathcal{K}_{\text{DG}}}},$
(V.22)

so the effective Schmidt number naturally arises and – aside from the benign
factor of $\sqrt{2}$ – identifies the weakly or strongly squeezed limit as
$|\mathring{\beta}|\ll 1$ or $|\mathring{\beta}|\gg 1$ respectively,
justifying the definition in section III. Since each $\overline{\chi}_{n}(t)$
has a width $\tau=\mathcal{T}_{c}^{\text{DG}}$, the effective Schmidt number
$\mathcal{K}_{\text{DG}}$ roughly identifies the number of Whittaker-Shannon
modes that are relevant along $r_{nn}$. Further, following the discussion
around Eq. (IV.19), for a very long pulse such that
$\mathcal{K}_{\text{DG}}\gg 1$ the squeezing parameter $|\beta|$ can be quite
large but $|\mathring{\beta}|$ remains finite.

In Fig. 14 we plot $|\overline{\gamma}(t_{1},t_{2})|^{2}$ and
$|\gamma(\omega_{1},\omega_{2})|^{2}$, reconstructed from using the Whittaker-
Shannon decomposition for the double-Gaussian joint amplitude of Fig. 2, where
we have chosen $\Omega=2\pi\mathcal{B}_{c}$, and plot $r_{nm}$. Comparing to
Fig. 2, we see very good agreement, and from the zoomed in plot of $r_{nm}$ in
the lower right corner it is clear that only a few neighbouring Whittaker-
Shannon modes in the $|n-m|$ direction are relevant.

The argument that even though $|\beta|$ can be quite large
$|\mathring{\beta}|$ remains finite holds true for _any_ joint amplitude,
because as long as it is square normalized (I.3) it will carry pre-factors on
its behavior in the two directions in the plane (of either $(t_{1},t_{2})$ or
$(\omega_{1},\omega_{2})$) over which it is defined. In Appendix C, we show
that given a joint temporal amplitude characterized by $\mathcal{T}_{p}$ and
$\mathcal{T}_{c}$, its maximum value $\overline{\gamma}_{\text{max}}$ is on
the order of

$\overline{\gamma}_{\text{max}}\sim\frac{1}{\sqrt{\mathcal{T}_{p}\mathcal{T}_{c}}}$
(V.23)

To apply the Whittaker-Shannon decomposition we set
$\Omega=2\pi\mathcal{B}_{c}$ so that $\tau=\mathcal{T}_{c}$. Then one
immediately finds

$\begin{split}|\mathring{\beta}|=\tau|\beta|\overline{\gamma}_{\text{max}}\sim\frac{|\beta|}{\sqrt{\mathcal{K}}}.\end{split}$
(V.24)

What remains to be shown is the relation between the Schmidt number $K$ and
the effective Schmidt number $\mathcal{K}$. In Appendix D, we show that for a
general joint temporal amplitude characterized by widths $\mathcal{T}_{p}$,
$\mathcal{T}_{c}$ indicated schematically in Fig. 1, and a given $\tau$ which
we set to be equal to $\mathcal{T}_{c}$, the inequality

$K\leq\frac{\mathcal{T}_{p}}{\tau}=\frac{\mathcal{T}_{p}}{\mathcal{T}_{c}}=\mathcal{K},$
(V.25)

generally holds; of course this is conditioned on $\tau$ being sufficiently
small that the Whittaker-Shannon decomposition can be accurately used. Thus by
applying the Whittaker-Shannon interpolation formula – which previously has
not been linked to squeezed light or the Schmidt decomposition – we are able
to show that the Schmidt number and effective Schmidt number are intimately
related. Further, the Whittaker-Shannon decomposition may be an important
stepping stone in formally linking the Schmidt number to other measures of the
correlation, like the time-bandwidth product and its generalizations, in the
same way that the Whittaker-Shannon interpolation formula is intimately linked
to the Shannon number for classical signals and the information content of
images [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 25, 26, 28, 29, 30, 31].

Finally, we are now in a position to explain why the effective Schmidt number
for the sinc-hat example is nearly equal to the exact Schmidt number. First we
emphasize that the sinc-hat is an idealization with a sharp cutoff in time
(frequency) and is constant along $t_{1}=t_{2}$ $(\omega_{1}=-\omega_{2})$.
This means that in the limit when $T_{p}\gg T_{c}$, to very good approximation
it is exactly bandlimited with $\Omega=\Omega_{c}$ ($\tau=T_{c}$).
Mathematically, the sharp cutoff in frequency results in an exact
diagonalization of $\overline{\gamma}(n\tau,m\tau)$, and since
$\overline{\gamma}(n\tau,n\tau)$ is identical for all $n$ where it is not zero
(see Eq. (IV.14)) the sum in Eq. (III.3) can be carried out exactly resulting
in $K\approx T_{p}/T_{c}$, as we found in section III.2; we then expect

$\frac{\mathcal{K}_{\mathrm{SH}}}{K_{\mathrm{SH}}}=1+\frac{T_{c}}{T_{p}}\to
1,$ (V.26)

for $T_{p}/T_{c}\to\infty$. Physically, the Schmidt amplitudes are near-
degenerate, so there is no unique set of Schmidt modes.

Figure 13: For the double-Gaussian, from left to right we plot the: joint
temporal intensity divided by its maximum value with the axes normalized by
$\mathcal{T}_{p}^{\mathrm{DG}}$; joint spectral intensity divided by its
maximum value with the axes normalized by $2\pi\mathcal{B}_{c}^{\mathrm{DG}}$;
amplitudes $r_{nm}$.

##  VI Employing the Whittaker-Shannon decomposition

In using the pseudo-Schmidt decomposition to “take apart” the sinc-hat joint
temporal amplitude, we showed that the expressions for the correlation
functions $\overline{G}^{(1)}(t)$ and $\overline{G}^{(2)}(t_{1},t_{2})$ could
be easily understood in terms of the properties of the supermodes involved in
the decomposition. In this section we look at the corresponding expressions
for the correlation functions when we “take apart” the joint amplitudes using
the Whittaker-Shannon decomposition instead. The results are more complicated,
but again the behavior of the correlation functions can be understood in terms
of the properties of the Whittaker-Shannon supermodes in an intuitive way. And
since the Whittaker-Shannon decomposition can be much more widely applied than
a pseudo-Schmidt decomposition, the results here are much more general.

To simplify the notation we write the squeezed state in Eq. (V.15) as
$\ket{\Psi}=\tilde{S}\ket{\text{vac}}$, where

$\tilde{S}=e^{\frac{1}{2}\sum\limits_{n,m}\beta_{nm}B_{n}^{\dagger}B_{m}^{\dagger}-\text{h.c.}},$
(VI.1)

is the squeezing operator. To calculate the correlation functions analogously
to what was done with the Schmidt and pseudo-Schmidt decompositions, we use
the inverse relation (Eq. (V.13)) and the transformation [32],

$\displaystyle\tilde{S}^{\dagger}B_{r}\tilde{S}=\mu_{rs}B_{s}+\nu_{rs}B_{s}^{\dagger},$
(VI.2a)
$\displaystyle\tilde{S}^{\dagger}B_{r}^{\dagger}\tilde{S}=\mu_{rs}^{*}B_{s}^{\dagger}+\nu_{rs}^{*}B_{s},$
(VI.2b)

where we adopt the convention that repeated indices are to be summed over, and

$\displaystyle\mu_{rs}=\delta_{rs}+\frac{1}{2!}\beta_{ra}\beta_{as}^{*}+\frac{1}{4!}\beta_{ra}\beta_{ab}^{*}\beta_{bc}\beta_{cs}^{*}+...,$
(VI.3a)
$\displaystyle\nu_{rs}=\beta_{rs}+\frac{1}{3!}\beta_{ra}\beta^{*}_{ab}\beta_{bs}+...\,.$
(VI.3b)

Note that from the symmetric property of $\beta_{nm}$ it follows that
$\mu_{rs}$ is Hermitian and $\nu_{rs}$ is symmetric.

For the Schmidt or pseudo-Schmidt decomposition the transformation used always
involves a single supermode (see Eq. (III.17)); for the general Whittaker-
Shannon decomposition the squeezing transformation is more complicated. The
structure of the squeezing transformation (VI.2), and the form of $\mu_{rs}$
and $\nu_{rs}$, motivates the use of matrix multiplication; $\beta_{nm}$ is
now treated as a complex square symmetric matrix which we denote by
$\bm{\beta}$ ($\bm{\beta}^{*}=\bm{\beta}^{\dagger}$). To implement the
squeezing transformation it is convenient to use the “left” polar
decomposition of $\bm{\beta}$ (valid for _any_ complex square matrix), given
by

$\bm{\beta}=\mathbf{UP},$ (VI.4)

where $\mathbf{P}=\sqrt{\bm{\beta^{\dagger}\beta}}$ is Hermitian and
$\mathbf{U}$ is unitary. Equivalently, if we set
$\mathbf{Q}=\mathbf{UPU}^{\dagger}$ then we have the “right” polar
decomposition $\bm{\beta}=\mathbf{QU}$; since $\bm{\beta}$ is symmetric,
$\mathbf{P}^{T}=\mathbf{Q}$ and $\mathbf{U}^{T}=\mathbf{U}$. Using the polar
decomposition and its properties, the matrices $\bm{\mu}$ and $\bm{\nu}$ are
given by [33]

$\displaystyle\bm{\mu}=\text{cosh}(\mathbf{UPU^{\dagger}})=\text{cosh}\mathbf{Q},$
(VI.5a)
$\displaystyle\bm{\nu}=\mathbf{U}\text{sinh}\mathbf{P}=(\text{sinh}\mathbf{Q})\mathbf{U}.$
(VI.5b)

The form of $\bm{\mu}$ and $\bm{\nu}$ guarantees that the transformation in
Eq. (VI.2) preserves the commutation relations of the $B_{n}$ operators, as
expected since the transformation is unitary.

In Appendix E we calculate the correlation functions (Eq. (II.4)); we find

$\begin{split}\overline{G}^{(1)}(t_{1},t_{2})&=\bm{\overline{\chi}}^{\dagger}(t_{1})(\text{sinh}^{2}\mathbf{P})\bm{\overline{\chi}}(t_{2}),\end{split}$
(VI.6)

and

$\begin{split}\overline{G}^{(2)}(t_{1},t_{2})&=|\bm{\overline{\chi}}^{T}(t_{1})\mathbf{U}(\text{sinh}\mathbf{P})(\text{cosh}\mathbf{P})\bm{\overline{\chi}}(t_{2})|^{2}\\\
&+\overline{G}^{(1)}(t_{1})\overline{G}^{(1)}(t_{2})+|\overline{G}^{(1)}(t_{1},t_{2})|^{2}\end{split}$
(VI.7)

where
$\bm{\overline{\chi}}(t)=(...,\overline{\chi}_{-1}(t),\overline{\chi}_{0}(t),\overline{\chi}_{1}(t),...)^{T}$
is the column vector formed from the set $\\{\overline{\chi}_{n}(t)\\}$ for a
given $t$.

Since the functions $\overline{\chi}_{n}(t)$ are similar to the pseudo-Schmidt
modes $\overline{\eta}_{n}(t)$, Eqs. VI.6 and VI.7 for the correlation
functions calculated using the Whittaker-Shannon decomposition are the
generalization of the pseudo-Schmidt results (IV.16, IV.17), valid for an
approximately bandlimited, but otherwise general joint amplitude. Both the
pseudo-Schmidt and Whittaker-Shannon mode functions are localized, so for
short time differences the structure of Eq. VI.7 reduces to that of the
pseudo-Schmidt decomposition, and we can again think of the correlation
function on a “mode-by-mode” basis; however, since this correspondence is only
approximate, we discuss it in Appendix F.

###  VI.1 Packet expansion

From Eq. (VI.6), we can immediately write the photon density as

$\begin{split}\overline{G}^{(1)}(t)&=\sum_{n}\Gamma_{n}^{2}|\overline{\rho}_{n}(t)|^{2},\end{split}$
(VI.8)

where $\overline{\rho}_{n}(t)$ is a normalized function set by

$\overline{\rho}_{n}(t)=\frac{1}{\Gamma_{n}}\sum_{m}(\text{sinh}\mathbf{P})_{nm}\overline{\chi}_{m}(t),$
(VI.9)

and

$\Gamma_{n}=\sqrt{(\text{sinh}^{2}\mathbf{P})_{nn}},$ (VI.10)

which is real. The expression (VI.8) for $\overline{G}^{(1)}(t)$ is the
generalization of Eq. (IV.16) in the pseudo-Schmidt decomposition, and clearly
has the same form; indeed, if we were to set $\beta_{nm}$ to be diagonal and
independent of $n$ for the $n$ for which it does not vanish, we would have
$\overline{\rho}_{n}(t)\to\overline{\chi}_{n}(t)$. Even more generally, the
expression (VI.8) mirrors the form of the expansion (III.20) of
$\overline{G}^{(1)}(t)$ in terms of Schmidt modes, with $\Gamma_{n}$ here
taking the role of $s_{n}$ there. But the $\overline{\rho}_{n}(t)$ cannot be
identified as a supermode; while the functions in the set
$\\{\overline{\chi}_{n}(t)\\}$ are mutually orthogonal, the functions in the
set $\\{\overline{\rho}_{n}(t)\\}$ are not, because in general
$\text{sinh}\mathbf{P}$ is not unitary. Nonetheless, the functions in the
latter set are generally localized compared to the duration of the pulse,
especially for weak squeezing. We refer to the $\overline{\rho}_{n}(t)$ simply
as “packets,” and to the expansion (VI.8) for $\overline{G}^{(1)}(t)$ as its
“packet expansion;” we will see packet expansions of other correlation
functions below.

The expected number of photons in the pulse is given by integrating
$\overline{G}^{(1)}(t)$ over all time; we find

$N_{\text{pulse}}=\text{Tr}(\text{sinh}^{2}\mathbf{P})=\sum_{n}\Gamma_{n}^{2},$
(VI.11)

where $\text{Tr}(\cdot)$ denotes the trace. This is reminiscent of the
corresponding expressions (III.22) and (IV.21) for the Schmidt and pseudo-
Schmidt expansions respectively. From Eq. (VI.11) it is clear that
$\Gamma_{n}^{2}$ is the number of photons in each packet, and summing over all
packets gives the total number of photons.

For the double-Gaussian and three values of $\beta$ chosen in section III.1,
we have $|\mathring{\beta}|\approx 0.014,0.7,$ and $1.4$ which is on the order
of the three values $|\beta|/\sqrt{\mathcal{K}}$, in agreement with the
discussion surrounding Eq. (V.24). In Fig. 14 we plot $G^{(1)}(t)$ calculated
using Eq. (VI.6), as well as the contributions given by Eq. (VI.8) for a few
values of $n$, together with the exact $\overline{G}^{(1)}(t)$ for the three
chosen values of $\beta$, which correspond to $N_{\text{pulse}}\approx
0.01,35,383$; we find excellent agreement between the exact and Whittaker-
Shannon decomposition. From Fig. 14 we see that each $\overline{\rho}_{n}(t)$
is clearly localized compared to the duration of the pulse, and so using the
packets we can “take apart” the squeezed light and provide a simple
description of the photon density (Eq. (VI.8)); this extends our understanding
from the pseudo-Schmidt decomposition valid for the sinc-hat joint amplitude
to more general joint amplitudes, such as the double-Gaussian, where a
Whittaker-Shannon decomposition is necessary.

Notice that as $|\beta|$ ($|\mathring{\beta}|$) increases so does the width of
each $\overline{\rho}_{n}(t)$. Referring back to Eq. (VI.9), this occurs
because elements of $\text{sinh}\mathbf{P}$ that are further off-diagonal
become more important as $|\beta|$ increases. And this is a consequence of the
fact that more powers of $\mathbf{P}$ become important in the expansion of
$\text{sinh}\mathbf{P}$ as $|\beta|$ increases, since $\mathbf{P}$ depends on
$\bm{\beta}$ (see Eq. (VI.4)). Thus the off-diagonality of
$\text{sinh}\mathbf{P}$ is extended beyond that of $\bm{\beta}$, and elements
of $\text{sinh}\mathbf{P}$ further from the diagonal become larger as
$|\beta|$ increases; see Fig. 15 for plots of $\text{sinh}\mathbf{P}$ with
increasing $|\beta|$, which demonstrates this effect.

Figure 14: For the double-Gaussian, from left to right we plot
$\overline{G}^{(1)}(t)/\Phi$ calculated using the Whittaker-Shannon
decomposition and a few contributions from different packets in Eq. (VI.8)
compared with the exact calculation (III.20), with the horizontal axis
normalized by $\mathcal{T}_{p}^{\mathrm{DG}}$ for $\beta=0.1,5,$ and $10$.
Figure 15: For the double-Gaussian, from left to right we plot
$(\text{sinh}\mathbf{P})_{nm}$ for $\beta=0.1,5,$ and $10$ with the horizontal
axis normalized by $\mathcal{T}_{p}^{\mathrm{DG}}$.

We can now construct the expansions for
$\overline{G}^{(2)}_{\text{coh}}(t_{1},t_{2})$ and
$\overline{G}^{(2)}_{\text{incoh}}(t_{1},t_{2})$. For the first of these,
comparing the expression (VI.7) for the full $\overline{G}^{(2)}(t_{1},t_{2})$
with our earlier general expression (III.32) for
$\overline{G}^{(2)}_{\text{incoh}}(t_{1},t_{2})$, we can identify

$\begin{split}\overline{G}^{(2)}_{\text{coh}}(t_{1},t_{2})&=|\bm{\overline{\chi}}^{T}(t_{1})\mathbf{U}(\text{sinh}\mathbf{P})(\text{cosh}\mathbf{P})\bm{\overline{\chi}}(t_{2})|^{2}\\\
&=\left|\sum_{n,m}(\mathbf{U}(\text{sinh}\mathbf{P})(\text{cosh}\mathbf{P}))_{nm}\overline{\chi}_{n}(t_{1})\overline{\chi}_{m}(t_{2})\right|^{2}\end{split}.$
(VI.12)

This can be compared with the corresponding expressions (III.21) and (IV.17)
for the Schmidt and pseudo-Schmidt decompositions respectively. The appearance
of terms with $m\neq n$ here, as opposed to the single summation that appears
in the Schmidt and pseudo-Schmidt expansion, is expected given that the
squeezed state written using the Whittaker-Shannon decomposition involves a
double sum (VI.1)

We show in Appendix G that the expression (VI.12) for
$\overline{G}^{(2)}_{\text{coh}}(t_{1},t_{2})$ can also be written using the
set of packet functions $\\{\overline{\rho}_{n}(t)\\}$,

$\overline{G}^{(2)}_{\text{coh}}(t_{1},t_{2})=\left|\sum_{n,m}\Gamma_{n}\Gamma_{m}(\mathbf{U}\text{coth}\mathbf{P})_{mn}\overline{\rho}_{m}(t_{1})\overline{\rho}_{n}(t_{2})\right|^{2},$
(VI.13)

but this is not a convenient expression to use in practice: For if
$\bm{\beta}$ is close to diagonal some functions of $\bm{\beta}$, such as
$\text{tanh}\mathbf{P}$, will be as well, but not $\text{coth}\mathbf{P}$.
Further, the weakly squeezed limit is not directly apparent from the form of
Eq. (VI.13), so it seems preferable to write
$\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})$ using
$\\{\overline{\chi}_{n}(t)\\}$ instead of $\\{\overline{\rho}_{n}(t)\\}$.
Perhaps this is not surprising, for earlier we found that when $|\beta|\ll 1$,
$\overline{G}^{(2)}(t_{1},t_{2})\to|\beta|^{2}|\overline{\gamma}(t_{1},t_{2})|^{2}$
(III.14), which directly depends on the correlations contained in the joint
temporal amplitude; in same limit, using the Whittaker–Shannon decomposition,
we expect it to depend on the analogous quantity $\beta_{nm}$. So unlike
$G^{(1)}(t)$, where the photon density at a particular time involves
contributions from all possible pairs and is written in terms of
$\\{\overline{\rho}_{n}(t)\\}$, $\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})$
– at least in the weakly squeezed limit – should directly depend on the
temporal correlations, and so it is more suitable to write
$\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})$ in terms of
$\\{\overline{\chi}_{n}(t)\\}$, as is done in Eq. (VI.12).

Turning finally to the general expression (III.32) for
$\overline{G}^{(2)}_{\text{incoh}}(t_{1},t_{2})$, using the definition of
$\overline{\rho}_{n}(t)$ and $\Gamma_{n}$ (VI.9, VI.10) we can write a packet
expansion for the incoherent contribution as

$\begin{split}\overline{G}^{(2)}_{\text{incoh}}(t_{1},t_{2})\\\
&=\frac{1}{2}\sum_{n,m}\left|\Gamma_{n}\Gamma_{m}(\overline{\rho}_{n}(t_{1})\overline{\rho}_{m}(t_{2})+\overline{\rho}_{n}(t_{2})\overline{\rho}_{m}(t_{1}))\right|^{2},\end{split}$
(VI.14)

which has the same form as the Schmidt (III.21) and pseudo-Schmidt (IV.17)
decompositions, but is in terms of the set of packets
$\\{\overline{\rho}_{n}(t)\\}$.

In Fig. 16 we plot $G^{(2)}(t_{1},t_{2})$ calculated using Eq. (VI.7) for the
double-Gaussian joint amplitude, as well as the coherent and incoherent
contributions; comparing to Fig. 4 for the exact calculation we find excellent
agreement between the two. Although the Whittaker-Shannon decomposition does
not allow the simple factorization of the ket into product kets associated
with each Schmidt or pseudo-Schmidt mode, correlation functions can still be
evaluated. And, since the off-diagonal elements of the squeezing matrix
$\bm{\beta}$ typically drop off quickly away from the diagonal, the
correlation functions can be written in a form involving either the set of
functions $\\{\overline{\chi}_{n}(t)\\}$ or the set of functions
$\\{\overline{\rho}_{n}(t)\\}$; all these functions are localized compared to
the Schmidt modes.

Figure 16: For the double-Gaussian and using Whittaker–Shannon decompositions,
from left to right we plot $\overline{G}^{(2)}(t_{1},t_{2})/\Phi^{2}$ (top)
and the coherent and incoherent contribution to
$\overline{G}^{(2)}(t/2,-t/2)/\Phi^{2}$ (bottom) with the axes normalized by
$\mathcal{T}_{p}^{\mathrm{DG}}$ for $\beta=0.1,5,$ and $10$.

###  VI.2 Correlation functions in the weakly squeezed limit

In this section we identify approximate expressions for the correlation
functions valid in the weakly squeezed limit when $|\mathring{\beta}|\ll 1$.
The correlation functions for the Whittaker–Shannon decomposition involve the
matrix $\mathbf{P}$, which using the form of Eq. (V.20) is given by
$\mathbf{P}=|\mathring{\beta}|^{2}\sqrt{\mathbf{r^{\dagger}r}}$. Taking the
weakly squeezed limit we approximate

$\begin{split}&\text{sinh}\mathbf{P}\to\mathbf{P}\\\
&\text{cosh}\mathbf{P}\to\mathbf{1},\end{split}$ (VI.15)

where $\mathbf{1}$ is the identity matrix. Then
$\Gamma_{n}\to\sqrt{(\mathbf{P}^{2})_{nn}}$,

$\overline{\rho}_{n}(t)\to\frac{1}{\Gamma_{n}}\sum_{m}P_{nm}\overline{\chi}_{m}(t),\\\
$ (VI.16)

and for $r_{nm}$ nonzero only for $|n-m|$ less than a small integer, the set
{$\overline{\rho}_{n}(t)\\}$ will be localized. The expression (VI.8) then
gives

$\begin{split}\overline{G}^{(1)}(t)&\to\sum_{n,m}|P_{nm}\overline{\chi}_{m}(t)|^{2}.\end{split}$
(VI.17)

Turning to $\overline{G}^{(2)}(t_{1},t_{2})$, for the coherent contribution
$\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})$, the general expression
(VI.12), using the result

$\mathbf{U}(\text{sinh}\mathbf{P})(\text{cosh}\mathbf{P})\to\mathbf{U}\mathbf{P}=\bm{\beta},$
(VI.18)

to find

$\begin{split}\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})&\to|\bm{\overline{\chi}}^{T}(t_{1})\bm{\beta}\bm{\overline{\chi}}(t_{2})|^{2}\\\
&=\left|\sum_{n,m}\overline{\chi}_{n}(t_{1})\beta_{nm}\overline{\chi}_{m}(t_{2})\right|^{2},\end{split}$
(VI.19)

which clearly shows that the resulting photon statistics depends on
$\beta_{nm}$. Finally, in this limit the expression (VI.14) for
$\overline{G}^{(2)}_{\text{incoh}}(t_{1},t_{2})$ gives

$\begin{split}\overline{G}^{(2)}_{\text{incoh}}(t_{1},t_{2})\to\\\
&\frac{1}{2}\sum_{n,m}\left|\sum_{u,v}P_{nu}P_{mv}(\overline{\chi}_{u}(t_{1})\overline{\chi}_{v}(t_{2})+\overline{\chi}_{u}(t_{2})\overline{\chi}_{v}(t_{1}))\right|^{2}.\end{split}$
(VI.20)

##  VII Local states and correlation functions

The general equations we derived for $\overline{G}^{(1)}(t)$ (VI.8) and
$\overline{G}^{(2)}(t_{1},t_{2})$ (VI.7), and their weakly squeezed
approximations, are valid for any times $t,t_{1},t_{2}$. In the discussion
surrounding Fig. 14 for the photon density, it was noted that since each
packet is localized we only need a few to properly represent the photon
density at any particular time. This suggests that for some calculations,
including some more general than the correlation functions considered above,
we can rely on an approximate form of the ket itself.

Suppose we are interested in features of the state around a small
neighbourhood centered at the time $t_{\textup{I}}$. Now for a correlated
joint temporal amplitude $\beta_{nm}$ will typically only be nonzero for
$|n-m|$ ranging up to a small integer, so for times around a small
neighbourhood of $t_{\textup{I}}$ only a few Whittaker-Shannon modes centered
around $n_{\textup{I}}=[t_{\textup{I}}/\tau]$ will be relevant; here we use
$[\cdot]$ to denote the closest integer. We identify the range of $m$ around
$n_{\textup{I}}$ for which $\beta_{n_{\textup{I}}m}$ will be non-negligible by
the odd integer $d$, assuming that $\beta_{n_{\textup{I}}m}$ is sufficiently
small for $|n_{\textup{I}}-m|>(d-1)/2$ that for those values of $m$ it can be
neglected.

We then split the matrix $\bm{\beta}$ into two contributions

$\bm{\beta}=\mathbf{R}^{\textup{I}}+\mathbf{K},$ (VII.1)

where $\mathbf{R}^{\textup{I}}$ is a symmetric matrix with nonzero entries
centered at $\mathrm{R}^{\textup{I}}_{n_{\textup{I}}n_{\textup{I}}}$. It
contains the elements of $\bm{\beta}$ as the row and column indices range over
$(d-1)/2$ in all directions from the center at
(${n_{\textup{I}}},{n_{\textup{I}}}$), and all its other elements are set to
zero. The matrix $\mathbf{R}^{\textup{I}}$ is shown schematically in Fig. 17
with the nonzero contributions existing inside the red “box” containing
$d^{2}$ elements; all the other elements of $\bm{\beta}$ are contained in
$\mathbf{K}$. For times of interest we assume that $d$ is chosen large enough
so that significant contributions to the quantities of interest, such as
correlation functions involving times near $t_{\textup{I}}$, only involve the
elements of $\bm{\beta}$ contained in $\mathbf{R}^{\textup{I}}$.

Figure 17: Schematic of the matrix $\mathbf{R}^{\textup{I}}$ which has nonzero
entries centered at $\beta_{n_{\textup{I}}n_{\textup{I}}}$ with a width $d$
and zeros everywhere else. The matrix
$\mathbf{K}=\bm{\beta}-\mathbf{R}^{\textup{I}}$ and consists of every nonzero
element that we set to zero in $\mathbf{R}^{\textup{I}}$. Figure 18: For the
double-Gaussian, from left to right we plot $\overline{G}^{(1)}(t)/\Phi$
calculated using the Whittaker-Shannon decomposition with the full
$\beta_{nm}$ compared to the approximate calculation using $R^{\textup{I}}$
near $t_{\textup{I}}=0$ for $d=7,9$ and $11$, with the horizontal axis
normalized by $\mathcal{T}_{p}^{\mathrm{DG}}$ for $\beta=0.1,5,$ and $10$.
Figure 19: Schematic of the matrix $\bm{\beta}$ partitioned into a set of non-
overlapping matrices $\mathbf{R}^{J}$, each with nonzero values centered at
$\beta_{n_{J}n_{J}}$ of size $d_{J}$ denoted by the red squares. The matrix
$\mathbf{K}=\bm{\beta}-\mathbf{R}^{\textup{I}}-\mathbf{R}^{\textup{II}}-...$,
and consists of every other non zero element contained in $\bm{\beta}$.

The squeezing operator (VI.1) is then approximated by

$\begin{split}\tilde{S}&=e^{\frac{1}{2}\sum\limits_{n,m}(R^{\textup{I}}_{nm}+K_{nm})B_{n}^{\dagger}B_{m}^{\dagger}-\text{h.c.}}\\\
&\approx
e^{\frac{1}{2}\sum\limits_{n,m}R^{\textup{I}}_{nm}B_{n}^{\dagger}B_{m}^{\dagger}-\text{h.c.}},\end{split}$
(VII.2)

where in the second line we dropped all the terms from $\mathbf{K}$. The
dimension of $\mathbf{R}^{\textup{I}}$ can be quite large if the pulse is
long, but since most of its entries are zero and the only nonzero entries have
a width given by $d$, we use a “prime” symbol below to indicate that we can
restrict the sum to be only over the nonzero elements of
$\mathbf{R}^{\textup{I}}$. In practice, this can drastically increase the
efficiency of numerical computations involving only a limited region of time.
The state is then taken to be
$\ket{\Psi_{\textup{I}}}=\tilde{S}_{\textup{I}}\ket{\text{vac}}_{\textup{I}}$,
where $\ket{\text{vac}}_{\textup{I}}$ is the vacuum state corresponding to the
modes associated with the nonzero elements in $\mathbf{R}^{\textup{I}}$, and
we set

$\tilde{S}_{\textup{I}}=e^{\frac{1}{2}\sum^{\prime}\limits_{n,m}R^{\textup{I}}_{nm}B_{n}^{\dagger}B_{m}^{\dagger}-\text{h.c.}}.$
(VII.3)

From the approximate state $\ket{\Psi_{\textup{I}}}$ we can apply the same
steps as above to calculate the correlation function, but instead of using the
full matrix $\bm{\beta}$ for the transformation and polar decomposition we use
the reduced matrix $\mathbf{R}^{\textup{I}}$. Then for times of interest we
have (see Eq. (V.13))

$\begin{split}\overline{a}(\overline{t}_{\textup{I}})&=\sum_{n}\overline{\chi}_{n}(\overline{t}_{\textup{I}})B_{n}\\\
&\approx\sum_{n}^{\prime}\overline{\chi}_{n}(\overline{t}_{\textup{I}})B_{n},\end{split}$
(VII.4)

where again we restrict the sum to be over the modes associated with the
nonzero elements of $\mathbf{R}^{\textup{I}}$, as in the approximate state
$\ket{\Psi_{\textup{I}}}$. The equations for the correlation functions (VI.6,
VI.7) can be applied with $\bm{\beta}$ replaced by $\mathbf{R}^{\textup{I}}$
for times near $t_{\textup{I}}$, with the appropriate restriction of the sums.

In Fig. 18 we plot a zoomed-in version of $\overline{G}^{(1)}(t)$ calculated
from the full state $\ket{\Psi}$ and the approximate state
$\ket{\Psi_{\textup{I}}}$ around the time $t_{\textup{I}}=0$ for the double-
Gaussian, the three values of $\beta$ ($\mathring{\beta}$), and $d=7,9$ and
$11$. We see that the choice of $d$ generally determines over how wide a
neighbourhood around $t_{\textup{I}}$ the contributions from the full state
are well approximated by the contributions from $\ket{\psi}_{\textup{I}}$. For
$\beta=0.1$ ($\mathring{\beta}=0.014$), which is well within the weakly
squeezed limit, only three Whittaker-Shannon modes on each side of the center
mode are relevant to accurately determine the photon density at $t=0$; this is
much less then the total number of modes along $r_{nn}$ which is set by
$\mathcal{K}$, and for this example is $\mathcal{K}_{\mathrm{DG}}=100$. Thus,
the state $\ket{\Psi_{\textup{I}}}$ provides a “local” description of the
photon density around $t=t_{\textup{I}}$. With the state
$\ket{\Psi_{\textup{I}}}$ one can also calculate $G^{(2)}(t_{1},t_{2})$ as
long as both times are near the time $t_{\textup{I}}$. We find similar
agreement and trends with $\beta$ and $d$ as for the photon density.

As $|\beta|$ ($|\mathring{\beta}|$) increases more neighbouring Whittaker-
Shannon modes are required to accurately reproduce the photon statistics at a
given time and we need to increase the size of the nonzero box of elements in
$\mathbf{R}^{\textup{I}}$. For as $|\beta|$ increases there is a larger
amplitude for photons to be described by different modes spread further apart
from each other; see the discussion in the paragraph after Eq. (VI.11). In the
same way that $(\text{sinh}\mathbf{P})_{nm}$ spreads in the $|n-m|$ direction
as $|\beta|$ increases (see Fig. 15) we need to choose a larger box to capture
all the possible contributions near a given $t_{\textup{I}}$. We return to
this point below.

Suppose now we are interested in the properties of the state associated with
two or more times $t_{\textup{I}},t_{\textup{II}},t_{\textup{III}},...$,
“sufficiently far apart” from one another. Then we can split $\bm{\beta}$ into
a _set_ of contributions given by

$\bm{\beta}=\mathbf{R}^{\textup{I}}+\mathbf{R}^{\textup{II}}+\mathbf{R}^{\textup{III}}+...+\mathbf{K},$
(VII.5)

with $\mathbf{R}^{\textup{I}}$ associated with $t_{\textup{I}}$,
$\mathbf{R}^{\textup{II}}$ associated with $t_{\textup{II}}$, etc., and where
by “sufficiently far apart” we mean that the corresponding boxes of sizes
$d_{\textup{I}},d_{\textup{II}}$, etc. associated with the regions of
$\mathbf{R}^{\textup{I}},\mathbf{R}^{\textup{II}}$, etc. that contain nonzero
elements do not overlap; this is shown schematically in Fig. 19, where we
indicate the regions of $\mathbf{R}^{\textup{I}},\mathbf{R}^{\textup{II}}$,
etc. that contain nonzero elements by I,II, etc. Again, the matrix
$\mathbf{K}$ contains the remaining contributions to $\bm{\beta}$ not in any
of the nonzero regions of the $\mathbf{R}^{J}$ matrices,
$J=\textup{I},\textup{II}$, etc. Then since each $\mathbf{R}^{J}$ has nonzero
elements only in the region where the others do not, each matrix in
$\\{\mathbf{R}^{J}\\}$ commutes with the rest and each contribution to the
squeezing operator can be split apart. So for the times of interest the state
is given by

$\ket{\Psi}\approx\bigotimes_{J}\ket{\Psi_{J}}=\bigotimes_{J}\tilde{S}_{J}\ket{\text{vac}}_{J},$
(VII.6)

where $\ket{\Psi_{J}}$, $\tilde{S}_{J}$, and $\ket{\text{vac}}_{J}$ are the
obvious generalization of $\ket{\Psi_{\textup{I}}}$, $\tilde{S}_{\textup{I}}$
and $\ket{\text{vac}}_{\textup{I}}$. Using equation (VI.11) for the average
photon number we similarly calculate that each time region $t_{J}$ has

$N_{J}=\text{Tr}(\text{sinh}^{2}\mathbf{P}^{J})$ (VII.7)

photons, where $\mathbf{P}^{J}$ is the matrix calculated from doing a polar
decomposition of the corresponding $\mathbf{R}^{J}$.

To calculate the correlation functions in the neighbourhood of a time $t_{J}$
we again use Eqs. VI.6 and VI.7 with the replacement of $\beta_{nm}$ with the
appropriate $\mathbf{R}^{J}$ as discussed above. If instead we want to
calculate $\overline{G}^{(2)}(t_{J},t_{J^{\prime}})$ for $J\neq J^{\prime}$
then the corresponding operators $a(t_{J})$ and $a^{\dagger}(t_{J^{\prime}})$
in Eq. (II.4) commute and the resulting second-order correlation function is

$\overline{G}^{(2)}(t_{J},t_{J^{\prime}})=\overline{G}^{(1)}(t_{J})\overline{G}^{(1)}(t_{J^{\prime}}),$
(VII.8)

where the photon densities evaluated at the times $t_{J},t_{J^{\prime}}$ are
evaluated using the respective contributions from $\ket{\Psi_{J}}$ and
$\ket{\Psi_{J^{\prime}}}$.

###  VII.1 Disentangling the squeezing operator

While the approximation of the state $\ket{\Psi}$ into the set of states
{$\ket{\Psi_{J}}$}, takes apart the squeezed state and provides a local
calculation of the correlation functions, it does not really give us intuition
of the state itself. To gain insight into that, we make use of the general
“disentangling formula” [33] applied to each squeezing operator
$\tilde{S}_{J}$,

$\tilde{S}_{J}=|\mathbf{W}^{J}|^{\frac{1}{2}}e^{\frac{1}{2}\sum\limits_{n,m}^{\prime}T_{nm}^{J}B_{n}^{\dagger}B_{m}^{\dagger}}e^{\sum\limits_{n,m}^{\prime}L_{nm}^{J}B_{n}^{\dagger}B_{m}}e^{-\frac{1}{2}\sum\limits_{n,m}^{\prime}T^{J,*}_{nm}B_{n}B_{m}},$
(VII.9)

where

$\displaystyle|\mathbf{W}^{J}|^{\frac{1}{2}}=\sqrt{\text{det}(\text{sech}(\mathbf{Q}^{J}))}$
(VII.10a)
$\displaystyle\mathbf{T}^{J}=(\mathbf{T}^{J})^{T}=\text{tanh}(\mathbf{Q}^{J})\mathbf{U}^{J}=\mathbf{U}^{J}\text{tanh}(\mathbf{P}^{J})$
(VII.10b)
$\displaystyle\mathbf{L}^{J}=\text{ln}(\text{sech}(\mathbf{Q}^{J})),$
(VII.10c)

and with $\mathbf{Q}^{J}$, $\mathbf{U}^{J}$ and $\mathbf{P}^{J}$ the same as
before (VI.4), but calculated from the reduced matrix $\mathbf{R}^{J}$. Then
acting the disentangled squeezing operator on the vacuum state we have for
each $\ket{\Psi_{J}}$

$\begin{split}\ket{\Psi_{J}}&=|\mathbf{W}^{J}|^{\frac{1}{2}}e^{\frac{1}{2}\sum\limits_{n,m}^{\prime}T^{J}_{nm}B_{n}^{\dagger}B_{m}^{\dagger}}\ket{\text{vac}}_{J}\\\
&\equiv\overline{S}_{J}\ket{\text{vac}}_{J},\end{split}$ (VII.11)

where $\overline{S}_{J}$ is the disentangled squeezing operator after acting
on the vacuum state. We point out that one could also apply the disentangling
formula to the whole state valid at all times, but this is not very
illuminating because one can already calculate the full correlation functions;
further, had we first applied the disentangling formula and then reduced the
state by getting rid of the terms that are negligible, the resulting state
would not be normalized, whereas $\ket{\Psi_{J}}$, as given in Eq. (VII.11)
always is.

Figure 20: For the double-Gaussian joint amplitude, we plot
$\overline{G}^{(1)}(t)/\Phi$ calculated using the Schmidt decomposition (top)
and Whittaker-Shannon decomposition (bottom) as well as a few contributions
from each calculation, for $\beta=150$ with the horizontal axis normalized by
$\mathcal{T}_{p}^{\mathrm{DG}}$. Figure 21: For the double-Gaussian joint
amplitude, we plot $\overline{G}^{(2)}(t_{1},t_{2})/\Phi^{2}$ (top), and the
coherent and incoherent contribution to
$\overline{G}^{(2)}(t/2,-t/2)/\Phi^{2}$ (bottom) calculated using the Schmidt
decomposition, for $\beta=150$ with the horizontal and vertical axis
normalized by $\mathcal{T}_{p}^{\mathrm{DG}}$.

Unlike the squeezing operator in its “entangled form” (VI.1), all the
operators in the exponent of Eq. (VII.11) are creation operators which always
commute so we can equivalently write the state for each $J$ as

$\ket{\Psi_{J}}=|\mathbf{W}^{J}|^{\frac{1}{2}}\bigotimes_{n,m}e^{\frac{1}{2}T^{J}_{nm}B_{n}^{\dagger}B_{m}^{\dagger}}\ket{\text{vac}}_{J}.$
(VII.12)

Unfortunately this form of the state is not as intuitive as the _single_
product state in the pseudo-Schmidt decomposition because for a given $n$ we
must include contributions from every other $m$ in the range of $T^{J}_{nm}$.

###  VII.2 The ket in the weakly squeezed limit

Instead, in situations where $|\beta|$ can be quite large but
$|\mathring{\beta}|$ is sufficiently small, as in a long pulse, we can take
advantage of the fact that it is $|\mathring{\beta}|$ that sets the magnitude
of the matrix $T^{J}_{nm}$. Given that the sum in Eq. (VII.11) is only over
the modes of interest, and not the whole joint amplitude, in the weakly
squeezed limit when $|\mathring{\beta}|\ll 1$ we can Taylor expand the
exponential so that

$\begin{split}\ket{\Psi_{J}}&\approx|\mathbf{W}^{J}|^{\frac{1}{2}}\left(\ket{\text{vac}}_{J}+\sqrt{\frac{N_{J}}{2}}\ket{\textup{II}}_{J}+...\right),\end{split}$
(VII.13)

where we have used
$(\mathbf{T}^{J})^{\dagger}\mathbf{T}^{J}\to(\mathbf{R}^{J})^{\dagger}\mathbf{R}$,
Eq. (VII.7) reduces to

$N_{J}\to\text{Tr}((\mathbf{R}^{J})^{\dagger}\mathbf{R}^{J}),$ (VII.14)

and we have introduced the normalized two-photon state for each $J$ by

$\ket{\textup{II}}_{J}=\frac{1}{\sqrt{2}}\sum_{nm}\frac{T_{nm}^{J}}{\sqrt{N_{J}}}B_{n}^{\dagger}B_{m}^{\dagger}\ket{\text{vac}}_{J}.$
(VII.15)

The state $\ket{\Psi_{J}}$ has $N_{J}\ll 1$ photons and the prefactor
$|\mathbf{W}^{J}|$ in the weakly squeezed limit is

$|\mathbf{W}^{J}|\approx 1-\frac{N_{J}}{2},$ (VII.16)

so the state remains normalized to order $N_{J}^{2}$, as expected. The two-
photon state $\ket{\textup{II}}_{J}$ is a superposition of all the ways in
which pairs of photons can be associated with the same Whittaker-Shannon
supermode, or different supermodes, within a neighbourhood of the time
$t_{J}$; one can easily extend the state in Eq. (VII.13) to higher order in
which two pairs, three pairs, etc. are considered. For a set of times,
$\\{t_{J}\\}$, in the weakly squeezed limit the full state in equation (V.15)
can be expanded as

$\ket{\Psi}\approx\bigotimes_{J}|\mathbf{W}^{J}|^{\frac{1}{2}}\left(\ket{\text{vac}}_{J}+\sqrt{\frac{N_{J}}{2}}\ket{\textup{II}}_{J}+...\right),$
(VII.17)

providing a localized description of squeezing light for correlated but
otherwise arbitrary joint amplitudes. In the long pulse limit, despite the
fact that $|\beta|\to\infty$, $|\mathring{\beta}|$ remains finite and we can
describe the light in the weakly squeezed limit as being composed of
approximately two-photons within a neighbourhood around each time $t_{J}$.

##  VIII The strongly squeezed limit

In this section we consider the strongly squeezed limit, where
$|\beta|/\sqrt{\mathcal{K}}\gg 1$, or equivalently $|\mathring{\beta}|\gg 1$.
The results we derived in section V and VI are valid for any approximately
bandlimited joint amplitude and for any squeezing parameter $\beta$. However,
following the discussion around Fig. 3, for the double-Gaussian with
$\mathcal{K}_{\mathrm{DG}}=100$ we found that as $|\beta|$ increased fewer
Schmidt modes were required to calculate the correlation functions, although
many were required to calculate the joint amplitude.

To explore this further, consider the double-Gaussian joint amplitude (Fig. 2,
with $\mathcal{K}_{\mathrm{DG}}=100$) but for $\beta=150$, corresponding to
$|\beta|/\sqrt{\mathcal{K}_{\mathrm{DG}}}\approx 15$ or
$|\mathring{\beta}|\approx 21$, well in the strongly squeezed regime. In Fig.
20 we plot $\overline{G}^{(1)}(t)$ calculated using the Schmidt and Whittaker-
Shannon decompositions. Clearly the first Schmidt mode is sufficient to
produce an accurate $\overline{G}^{(1)}(t)$, despite the fact that all Schmidt
modes are required to correctly calculate the joint amplitude. This is because
the correlation functions depend on $s_{n}$ and $c_{n}$ (recall Eqs.
(III.20),(III.21)), but when $|\beta|$ is large these scale exponentially;
since the Schmidt modes drop off as $n$ increases, $s_{0}\gg s_{1}$ and the
sums in Eq. (III.20) for the correlation functions are well approximated by
the $n=0$ term. For the Whittaker-Shannon decomposition we still find very
good agreement with the Schmidt calculation. However, as in Fig. 14, we find
each packet is significantly broadened. In the strongly squeezed limit, the
amplitude that many photon pairs will contribute is large, and so the
contribution of photons corresponding to two Whittaker-Shannon modes for which
$|n-m|\gg 1$ is significant; since the ${\overline{\rho}_{n}}(t)$ include
contributions from all other $m$ for a given $n$, they are necessarily
broader. More mathematically, in the strongly squeezed limit many matrix
multiplications are involved in calculating $\text{sinh}\mathbf{P}$, and so
(following the discussion in section VI.1 and surrounding Fig. 15, 18,) the
width of each packet is significantly broadened.

Since packets can then significantly overlap with a number of their
neighbours, the local description of the photon statistics and resulting state
breaks down. This is not surprising, given that the second-order correlation
function, calculated using the Schmidt decomposition and plotted in Fig. 21,
is completely uncorrelated, a local description to identify the photon
correlations is not necessary. Note that here we do not include the plot of
$\overline{G}^{(2)}(t_{1},t_{2})$ calculated using the Whittaker-Shannon
decomposition, because it is essentially identical to Fig. 21.

Finally, we point out that in Fig. 21 the incoherent contribution is
approximately twice the coherent contribution. If we restrict the sum in Eq.
(III.20) and (III.21) to the first term we have

$\displaystyle\overline{G}^{(1)}(t)\to s_{0}^{2}|\overline{f}_{0}(t)|^{2},$
(VIII.1a) $\displaystyle N_{\text{pulse}}\to s_{0}^{2},$ (VIII.1b)
$\displaystyle\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})\to
N_{\text{pulse}}^{2}|\overline{f}_{0}(t_{1})|^{2}|\overline{f}_{0}(t_{2})|^{2},$
(VIII.1c) $\displaystyle\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2})\to
2N_{\text{pulse}}^{2}|\overline{f}_{0}(t_{1})|^{2}|\overline{f}_{0}(t_{2})|^{2},$
(VIII.1d)

where we have used the fact that for large $|\beta|$, $c_{n}\approx s_{n}$; it
is clear from these expressions that
$\overline{G}_{\text{incoh}}^{(2)}(t_{1},t_{2})\approx
2\overline{G}_{\text{coh}}^{(2)}(t_{1},t_{2})$. Then putting them together we
have

$\overline{G}^{(2)}(t_{1},t_{2})=3N_{\text{pulse}}^{2}|\overline{f}_{0}(t_{1})|^{2}|\overline{f}_{0}(t_{2})|^{2}.$
(VIII.2)

Since the calculation in Appendix F was done for a pseudo-Schmidt mode, it can
be applied to this instance of a single Schmidt mode, so we identify
$3N_{\text{pulse}}^{2}$ as the expectation value of the number of ways of
“picking” two photons in the large $|\beta|$ limit from the first Schmidt
mode.

Figure 22: From left to right we plot the: joint temporal intensity; joint
spectral intensity; and Schmidt amplitudes generated from a dual-pump
spontaneous four-wave mixing process. The two pump functions are centered at
the wavelengths $\overline{\lambda}_{\text{P1}}=1.556$ $\mu$m and
$\overline{\lambda}_{\text{P2}}=1.547$ $\mu$m and each have temporal FWHM of
$100$ ns and an energy of $10^{3}$ pJ. The generated photons are centered at
$\overline{\lambda}_{\text{S}}=1.552$ $\mu$m
($\overline{\omega}_{\text{S}}/2\pi=193.164$ THz) and have a bandwidth on the
order of a GHz. The ring resonator has quality factors
$Q_{\text{P1}}=1529378$, $Q_{\text{P2}}=3844257$, and $Q_{\text{S}}=2704405$
for the three modes and a nonlinear coupling $\Lambda=5$ THz [4].

##  IX A final example

In this section we consider a realistic joint amplitude generated from a dual-
pump spontaneous four-wave mixing process, in a ring resonator system, when
time ordering effects and self- and cross-phase modulation are included [4].
In Fig. 22 we plot the joint intensity and the Schmidt amplitudes. The joint
intensity has widths $\mathcal{T}_{p}^{\mathrm{R}}=200\text{ ns}$ and
$\mathcal{B}_{c}^{\text{R}}\approx 1.03\text{ GHz}$
($\mathcal{T}_{c}^{\text{R}}\approx 0.97\text{ ns}$) corresponding to an
effective Schmidt number $\mathcal{K}_{\text{R}}\approx 206$, and the joint
amplitude has a Schmidt number $K_{\text{R}}\approx 108$, where we use “R” to
identify this “ring” calculation. The squeezing parameter for the generation
is $\beta=3.72$, corresponding to $N_{\text{pulse}}\approx 11$ photons; other
system parameters, such as the center wavelengths and pump duration, are given
in the caption of Fig. 22. Comparing Fig. 22 with Fig. 23 for the joint
temporal intensity calculated using the Whittaker-Shannon decomposition, with
$\Omega/2\pi=\mathcal{B}_{c}^{\text{R}}$ ($\tau=\mathcal{T}_{c}^{\text{R}}$),
we see excellent agreement.

In Fig. 24 we plot $\overline{G}^{(1)}(t)$ for the Schmidt (top) and
Whittaker-Shannon (bottom) decompositions, as well as a few contributions to
each. We see a dramatic difference: At any particular time a huge number of
Schmidt modes are required to capture the overall photon density, while in our
“packet” decomposition only a small range of packets are required to describe
the behaviour at any time.

In Fig. (25) we plot $\overline{G}^{(2)}(t_{1},t_{2})$ (top) and the coherent
and incoherent contributions to $\overline{G}^{(2)}(t/2,-t/2)$ (bottom)
calculated using the Whittaker-Shannon decomposition; here we only plot the
Whittaker-Shannon results because the Schmidt results look the same. The
contributions to $\overline{G}^{(2)}(t_{1},t_{2})$ – strong peak near
$t_{1}=t_{2}$ and smaller broad background – match that of the double-Gaussian
example for the three values of $\beta$. This is unsurprising because,
although $\beta=3.72$ in this example, $\mathring{\beta}=0.28$, and so the
state is weakly squeezed.

Since the state is weakly squeezed, the formalism provided in section VII can
be directly applied, providing a localized description of the pulse of light
into a set of states with different time labels, each containing approximately
two photons.

Figure 23: Joint temporal intensity calculated using the Whittaker-Shannon
decomposition. Figure 24: Plot of $\overline{G}^{(1)}(t)$ calculated using the
Schmidt decomposition (top) and Whittaker-Shannon decomposition (bottom).
Figure 25: Plot of $G^{(2)}(t_{1},t_{2})$ (top) and the coherent and
incoherent contribution to $G^{(2)}(t/2,-t/2)$ calculated using the Whittaker-
Shannon decomposition.

##  X Conclusion

In this article we have developed a formalism to describe squeezed light with
a large spectral-temporal correlation. As opposed to the usual strategy of
employing the Schmidt decomposition, we feel it makes the physics more
apparent. We began by characterizing general joint amplitudes by their
timewidth $\mathcal{T}_{p}$ and bandwidth $\mathcal{B}_{c}$ (or equivalently
the coherence time $\mathcal{T}_{c}=1/\mathcal{B}_{c}$). Using the double-
Gaussian joint amplitude as an example, we calculated the correlation
functions using the Schmidt decomposition and found that for weak squeezing
the form of $\overline{G}^{(2)}(t_{1},t_{2})$ matches that of the joint
temporal amplitude and reaches its maximum value when
$|t_{1}-t_{2}|<\mathcal{T}_{c}$; however, the Schmidt modes themselves extend
on a much broader time scale given by $\mathcal{T}_{p}$. When calculating the
correlation functions using the Schmidt decomposition, we found that a large
amount of interference is present; we cannot associate a single Schmidt mode,
or even a few, with a particular time.

Next we considered another example, the sinc-hat joint amplitude, which
demonstrated that this behaviour of the double-Gaussian is not unique. And
although it is somewhat artificial, the sinc-hat joint amplitude is
interesting in that its Schmidt amplitudes are nearly degenerate, allowing us
to construct an approximate pseudo-Schmidt decomposition where the pseudo-
Schmidt modes are displaced, localized “sinc” functions. Using this
decomposition we could immediately identify contributions contained in
$\overline{G}^{(1)}(t)$ or $\overline{G}^{(2)}(t_{1},t_{2})$ at a particular
time as arising from a single pseudo-Schmidt mode, allowing us to “take apart”
the photon statistics, elucidating the physics. We also demonstrated that the
weakly squeezed limit corresponds to $|\beta|/\sqrt{\mathcal{N}}\ll 1$, where
$\mathcal{N}$ is the effective Schmidt number that identifies the number of
pseudo-Schmidt modes required for the decomposition. This is useful because in
the long pulse limit, where $|\beta|\gg 1$ and $N_{\text{pulse}}\gg 1$, the
quantity $|\beta|/\sqrt{\mathcal{N}}$ remains finite; despite there being a
large number of photons in the pulse, the number of photons in each pseudo-
Schmidt mode can be relatively small.

To consider more general joint amplitudes, where there is a range of Schmidt
amplitudes, we generalized the pseudo-Schmidt decomposition to any
approximately bandlimited joint amplitude by using the Whittaker-Shannon
interpolation formula. While the exponent in the squeezing operator then
involves a double sum instead of the usual single sum, we can nevertheless
define a packet expansion where each packet typically has a duration short
compared to that of the pulse, with the packets thus analogous to the pseudo-
Schmidt modes. In general, if the squeezing is weak to moderate, the
correlation functions at a particular time are associated with only a few
packets, allowing us to “take apart” the squeezed light and the resulting
photon statistics. Finally we showed that if one is only interested in some
finite time regions that form part of the pulse duration, an effective ket can
be written as a product of kets associated with those time regions. Then in
the weakly squeezed limit, which we can take to be set by the limit
$|\mathring{\beta}|\ll 1$, there will be on average only a few photons within
each time region, although in the long pulse limit the total number of photons
will be very large.

In extensions of this work, we will consider two- and higher-mode squeezing,
which is fairly straightforward, and we will apply this formalism to quantum
optics-based experiments such as coincidence-accidental-detection ratios and
SU(1,1) interferometry. And instead of describing squeezed light in the
spectral-temporal domain, one could expand this formalism to describe squeezed
light with a large correlation between the photon wavevector and the conjugate
position variables in two, or three dimensions. Another interesting focus is
the formulation of relationships between the Schmidt number, the effective
Schmidt number discussed here, and the Shannon number (or time-bandwidth
product) and its generalizations. Here there is the opportunity to apply the
vast mathematical formalism that has already been developed to describe the
information content of classical signals and images [8, 9, 10, 11, 12, 13, 14,
15, 16, 17, 25, 26, 28, 29, 30, 31].

With the continuous advancement of custom engineered nonlinear optical
systems, the generation of nonclassical “tri-photon” states is slowly becoming
a reality [34]. Such states and their “troint amplitudes” are the
generalization of squeezed states and the joint amplitudes discussed here. To
characterize these states, generalizations of the Schmidt decomposition
(singular-value-decomposition) need to be employed. One such example is the
canonical polyadic decomposition (or CP decomposition) [35]; however, an
orthogonal decomposition is not guaranteed to exist. Other generalizations of
the Schmidt decomposition exist, such as the Tucker decomposition [35], which
gains orthogonality but loses the single-sum behaviour of the Schmidt
decomposition. An alternative description can be provided using the Whittaker-
Shannon interpolation formula. In the same way we generalized it to two-
dimensional functions in this paper, it can be directly applied to any number
of dimensions in a very straightforward way. Thus the formalism applied here
can easily be generalized to situations where other decompositions are not
possible.

###### Acknowledgements.

We thank Xanadu Quantum Technologies for allowing access to their proprietary
code repository used in section IX, and Luke Helt for answering related
questions. We also thank Marco Liscidini and Nicolás Quesada for valuable
discussions. This work was supported by the Natural Sciences and Engineering
Research Council of Canada (NSERC). C. D. acknowledges an Ontario Graduate
Scholarship.

##  Appendix A Field operator

For a quasi-1D structure, if $z$ is the direction in which light is
propagating, the electric field operator in the Heisenberg picture takes the
form [4]

$\begin{split}\bm{E}(x,y,z;t)&=\sum_{n}\int_{-\infty}^{\infty}\frac{dk}{\sqrt{2\pi}}\>\bm{e}_{n}(k;x,y)b_{n}(k)e^{ikz}e^{-i\omega
t}+h.c.,\end{split}$ (A.1)

where the dependence of $\omega$ on $k$ is identifed by the dispersion
relation, the $\bm{e}_{n}(k;x,y)$ is a properly normalized field profile for a
transverse mode $n$ propagating with $k$ at frequency $\omega$, and $b_{n}(k)$
is the associated lowering operator,
$\left[b_{n}(k),b_{m}^{\dagger}(k^{\prime})\right]=\delta_{nm}\delta(k-k^{\prime}).$
We assume that only one transverse mode is of interest and drop the index $n$,
and that only $k>0$ are of interest. Taking $\omega=vk>0,$ where $v$ is the
group velocity, we identify modes by their frequency $\omega$, putting
$c(\omega)=v^{-1/2}b(k)$ so that
$\left[c(\omega),c^{\dagger}(\omega^{\prime})\right]=\delta(\omega-\omega^{\prime})$
holds, and we can then write

$\displaystyle\bm{E}(x,y,z;t)\rightarrow\int_{0}^{\infty}\frac{d\omega}{\sqrt{2\pi}}\frac{\bm{e}(\frac{\omega}{v};x,y)}{v^{1/2}}c(\omega)e^{i\omega
z/v}e^{-i\omega t}+h.c.$ (A.2)

Assuming that over the frequency range of interest $\bm{e}(\omega/v;x,y)$
varies little from its value $\bm{e}(\omega_{o}/v;x,y)$ at a center frequency
$\omega_{o}$, we can write

$\begin{split}&\bm{E}(x,y,z;t)\rightarrow\frac{\bm{e}(\frac{\omega_{o}}{v};x,y)}{v^{1/2}}\int_{0}^{\infty}\frac{d\omega}{\sqrt{2\pi}}c(\omega)e^{i\omega
z/v}e^{-i\omega t}+h.c.\\\
&=\frac{\bm{e}(\frac{\omega_{o}}{v};x,y)}{v^{1/2}}e^{i\omega_{o}z/v}e^{-i\omega_{o}t}\int_{0}^{\infty}\frac{d\omega}{\sqrt{2\pi}}c(\omega)e^{i(\omega-\omega_{o})z/v}e^{-i(\omega-\omega_{o})t}+h.c.\end{split}$
(A.3)

Then putting $a(\omega-\omega_{o})\equiv c(\omega),$ the commutation relations
(Eq. (I.2)) hold and for the range of frequencies much less than $\omega_{o}$
we have

$\begin{split}\bm{E}(x,y,z;t)&=\frac{\bm{e}(\frac{\omega_{o}}{v};x,y)}{v^{1/2}}e^{i\omega_{o}z/v}e^{-i\omega_{o}t}\int_{-\infty}^{\infty}\frac{d\omega}{\sqrt{2\pi}}a(\omega)e^{i\omega
z/v}e^{-i\omega t}+h.c.\end{split}$ (A.4)

This gives

$\displaystyle\bm{E}(x,y,0;t)=\frac{\bm{e}(\frac{\omega_{o}}{v};x,y)}{v^{1/2}}e^{-i\omega_{o}t}\overline{a}(t)+h.c.,$
(A.5)

where

$\displaystyle\overline{a}(t)=\int_{-\infty}^{\infty}\frac{d\omega}{\sqrt{2\pi}}a(\omega)e^{-i\omega
t}.$ (A.6)

Now since the Schrödinger operator for $\bm{E}(x,y,z)$ is just the Heisenberg
operator $\bm{E}(x,y,z;0),$ for that Schrödinger operator we have

$\displaystyle\bm{E}(x,y,z)=\frac{\bm{e}(\frac{\omega_{o}}{v};x,y)}{v^{1/2}}e^{i\omega_{o}z/v}\overline{a}(-\frac{z}{v})+h.c.$
(A.7)

That is, the operator $\overline{a}(t)$ is associated with the electric field
at time $t$ and $z=0$, and as well with the electric field at zero time and
position $z=-vt$, as would be expected because of the propagation with group
velocity $v$.

##  Appendix B Schematic of sinc-hat joint intensity

Consider the schematic sinc-hat joint intensity shown in Fig. 26, where we
drop the sinc “tails” along the anti-diagonal (diagonal) direction for the
joint temporal (spectral) intensity. Of course, because of the sinc tails the
joint intensities extend to infinity in either direction, but for
$T_{p}/T_{c}$ sufficiently large these contributions are small enough that the
behaviour in the anti-diagonal (diagonal) direction is effectively captured by
the width $T_{c}$ ($\Omega_{p}$) set at $t_{2}=0$ ($\Omega_{2}=0$).

In the schematic one can see that along the lines $t_{1}=t_{2}$
($\omega_{1}=-\omega_{2}$) the joint temporal (spectral) intensity ranges over
$-T_{p}/2\to T_{p}/2$ ($-\Omega_{c}/2\to\Omega_{c}/2$), however this is not
the _full_ horizontal extent of the joint intensities. For the joint temporal