# Not all Probability Density Functions are Tomograms

L.A. Markovich1,2,3,†, J. Urbanetz1 and V.I. Man’ko3,4
1 Institute-Lorentz, University Leiden, P.O. Box 9506,
2300 RA Leiden, The Netherlands
2 Institute for information transmission problems, Bol. Karetny per. 19,
Moscow 127051, Russia
3 Russian Quantum Center, Skolkovo, Moscow 121205, Russia
4 Lebedev Physical Institute, Russian Academy of Sciences,
Leninskii Prospect 53, Moscow 119991, Russia
∗Corresponding author e-mail<EMAIL_ADDRESS>
###### Abstract

This paper delves into the significance of the tomographic probability density
function (pdf) representation of quantum states, shedding light on the special
classes of pdfs that can be tomograms. Instead of using wave functions or
density operators on Hilbert spaces, tomograms, which are the true pdfs, are
used to completely describe the states of quantum systems. Unlike quasi-pdfs,
like the Wigner function, tomograms can be analysed using all the tools of
classical probability theory for pdf estimation, which can allow a better
quality of state reconstruction. This is particularly useful when dealing with
non-Gaussian states where the pdfs are multi-mode. The knowledge of the family
of distributions plays an important role in the application of both parametric
and non-parametric density estimation methods. We show that not all pdfs can
play the role of tomograms of quantum states and introduce the conditions that
must be fulfilled by pdfs to be ”quantum”.

##  1 Introduction

In the early 1900s, Schrödinger introduced the concept of a quantum system’s
state [1], characterised by a complex wave function. The notion of a density
matrix [2, 3, 4] or a density operator, acting in a Hilbert space, was
established shortly after to characterise the generic states of quantum
systems.

Nowadays, quantum state reconstruction is a critical task in quantum
information science and quantum technology. We know that in classical physics,
it is always possible to assign a unique probability distribution function
(pdf) to the outcomes of a physical measurement in continuous-variable
systems. Based on the measurement results, one can estimate the pdf using
different parametric and nonparametric methods [5] and thus fully describe the
system. However, in quantum mechanics, the Heisenberg uncertainty principle
states that certain pairs of non-commuting observables, such as position and
momentum, cannot be simultaneously precisely determined. Therefore, it is
challenging to reconstruct the quantum state from measurements, which leads to
the need for a more sophisticated representation of the quantum state.

Quasi-probability distributions (qpdf) are introduced to represent the
probabilities of outcomes of non-commuting observables, such as position and
momentum, in a way that respects the uncertainty principle. Qpdfs are not
actual probabilities in the classical sense. These distributions are often
complex-valued and can take negative values, which can make their
interpretation challenging. The Wigner function is one of the most well-known
qpdfs [6]. It provides a phase-space representation of a quantum state,
allowing the simultaneous description of position and momentum. The Wigner
function is real-valued and can take negative values. The Husimi Q-function is
another qpdf [7, 8] that is positive and can be thought of as a smoothed
version of the Wigner function [9]. The Glauber-Sudarshan P-function [10, 11]
is used to represent states in the context of quantum optics and can be used
for the description of coherent states. There exist other quasi-probabilities
defined as different Fourier-like transforms, like Kirkwood-Rihaczek (KR)
[12], Generalized KR [13], Margenau-Hill (MH) [14], Page [15], Choi-Williams
(CW) [16], sinc [17], and Cohen [18].

A symplectic tomogram [19, 20] is an alternative way to represent the quantum
state. It is a real-valued continuous probability distribution in the phase
space and is constructed through a convolution of the density operator with a
kernel that depends on the phase space coordinates. The resulting function
represents the probability of finding the quantum system in a particular
position and momentum configuration. The symplectic tomogram is related to the
qpdfs with the invertible integral transforms [21]. In fact, the Wigner
function can be derived from the symplectic tomogram through an inverse Radon
transform [22].

Since the symplectic tomogram is a pdf, one can use all the techniques known
in probability theory for pdf estimation. The parametric pdf estimation
methods strongly rely on the knowledge of the class of the pdf. For example,
the histogram method is widely used in practise [23, 24] because of its
simplicity. It divides the data range into intervals or bins and counts the
number of data points in each bin. The pdf is then approximated by normalizing
these counts. However, in some cases, a data sample does not resemble a common
probability distribution or cannot be easily made to fit the distribution.
This is often the case when the data has two peaks (bimodal distribution) or
many peaks (multi-modal distribution). In this case, the nonparametric methods
are used. These methods do not assume a specific parametric form for the pdf.
Hence, knowledge of the pdf family is essential for an optimal selection of
the pdf reconstruction method.

By definition, a state is called Gaussian if its Wigner function is a Gaussian
like function [25]. Since the Wigner function and the tomogram are connected
by the Radon transform, we can say that the state is Gaussian if its tomogram
is a Gaussian pdf. In [26] the quantum harmonic oscillator (HO) is studied,
showing that its tomograms are from the exponential family of distributions
[27]. In particular, the tomogram of the ground state is the Gaussian pdf. The
inverted oscillator [28] state tomogram is also from the exponential family.
The tomograms are known to be Gaussian pdfs for the squeezed vacuum states and
the thermal states [29]. However, a superposition of Gaussian states is
already a non-Gaussian mixture [30]. There has been a strong experimental
focus on two specific types of non-Gaussian states: cat states [31] and
Gottesman-Kitaev-Preskill (GKP) [32] states. In [33] the crystallized cat
states tomograms are presented. Being the superposition of cat states, their
states result in being non-Gaussian too since their tomogram is the sum of
Gaussian functions, which in general is not a Gaussian pdf.

Thus, non-Gaussian states are frequently encountered in practice. However, can
entirely different distributions describe a quantum state, for instance, from
the exponential or power-law family? It is evident that not all pdfs are
suitable for the role of a tomogram due to the constraints imposed on the
quantum state density operator, e.g., that it is nonnegative, Hermitian, and
has a unit trace. Thus, the question we address in this paper is: What
continuous distributions can describe quantum states? Knowing the answer, we
can solve the inverse problem: given a pdf playing the role of the tomogram,
find the new corresponding quantum state not known in the literature.

###  1.1 Contributions of this paper

In this manuscript, we introduce the characteristic function of the tomogram
pdf as a description of the quantum state. Being the Fourier transform of the
pdf, the characteristic function always exists and is bounded. We show that
the characteristic function of the pdf corresponding to the quantum state must
fulfil the conditions following from the density operator properties. If the
distribution does not satisfy these conditions, then the corresponding quantum
state does not exist. Thus, it is possible to assess the ability of known pdfs
to generate physical quantum states.

The trace of the density operator and the purity parameter are explicitly
written in terms of the characteristic function. Moreover, we show that the
trace of the product of two density operators can be written as an integral
from their characteristic functions corresponding to their tomogram pdfs.
Similarly the fidelity parameter is written. Thus, knowledge about the
characteristic function can be used for purity testing in the case of two
identical states or in the case of two different states for testing their
overlap. Also, the elements of the density matrix are written in terms of the
characteristic function, giving a novel way of parametrizating the density
matrix.

We investigate the exponential family of pdfs to find the conditions on its
parametrization to satisfy the constraints on the characteristic function. We
check some known pdfs from the exponential family: exponential, gamma, and
$\chi^{2}$ pdfs, showing that they never satisfy the latter conditions, thus
not being suitable for quantum state generation.

Finally, we study particular systems with continuous variables like
oscillators, the superposition of cat states, and different Gaussian states.
Their pdfs either belong to the exponential family of pdfs or are a mixture of
Gaussians. The corresponding characteristic functions are derived, showing
that they satisfy the conditions imposed on the characteristic function to
generate a physical quantum state. For a pseudoharmonic oscillator [34] that
is defined on the positive semiaxis, we deduce the tomogram for the first
time, showing it is from the exponential family too.

###  1.2 Organization of the paper

The paper is organised as follows. In Sec. 2 we recall the notion of the
symplectic tomogram, giving its connection to the density matrix of the
quantum state. In Sec. 3 the characteristic function of the tomogram is
introduced. The requirements that must be satisfied to characterize a quantum
state are outlined. In Sec. 4 the general exponential family of pdfs is
considered. We study the examples of continuous pdfs known in classical
probability theory. The quantum harmonic and pseudoharmonic oscillators are
studied in detail in Sec. 5. The superposition of cat states is analysed in
Sec. 6.

##  2 Symplectic Tomogram

Let us consider a quantum state in an infinite-dimensional Hilbert space
$\mathcal{H}$ associated with a positive Hermitian operator $\hat{\rho}$,
called the density matrix. The kernel of the density matrix in coordinate
representation is $\rho(x,x^{\prime})=\bra{x}\rho\ket{x^{\prime}}$, where
$\ket{x}$ is an eigensate of the position operator $\hat{q}$. The density
matrix operator is hermitian
($\rho^{\star}(x,x^{\prime})=\rho(x^{\prime},x)$). Its diagonal elements are
nonnegative ($\rho(x,x)\geq 0$) and its trace is equal to one
($\int\rho(x,x)dx=1$).

According to the general scheme [19], in the case of continuous variables, the
mapping of the density matrix to the family of pdfs, depending on two real
parameters $\mu$ and $\nu$, is given by the following relation:

$\displaystyle\mathcal{W}(X|\mu,\nu)=\text{Tr}{(\hat{\rho}\delta\left(\hat{1}X-\mu\hat{q}-\nu\hat{p}\right))}=\bra{X;\mu,\nu}\hat{\rho}\ket{X;\mu,\nu}.$
(1)

The $\delta(\cdot)$ is the Dirac delta-function. Here $\hat{q}$ and $\hat{p}$
are the position and momentum operators:

$\displaystyle\hat{q}\ket{X;\mu,\nu}=X\ket{X;\mu,\nu},$ (2)

and $|X;\mu,\nu\rangle$ is an eigenvector of the hermitian operator

$\displaystyle\hat{X}(\mu,\nu)=\mu\hat{q}+\nu\hat{p},\quad\hbar=1.$ (3)

The latter is a canonical transform of $\hat{q}$ and $\hat{p}$. Formally, this
quantity is a coordinate, measured in scaled and rotated reference frame in
the phase space. Its pdf (1) is called the symplectic tomogram. As a pdf the
tomogram is nonnegative and normalised:

$\displaystyle\int_{-\infty}^{\infty}\mathcal{W}(X|\mu,\nu)dX=1.$ (4)

For the symplectic tomogram, the inverse quantum Radon transform is the
following:

$\displaystyle\hat{\rho}=\frac{1}{2\pi}\int\mathcal{W}(X|\mu,\nu)\exp{(i(X\hat{1}-\mu\hat{q}-\nu\hat{p}))}dXd\mu
d\nu,$ (5)

defining the density matrix operator by the corresponding pdf. For the density
operator $\hat{\rho}=\ket{\psi}\bra{\psi}$ of the pure state $\ket{\psi}$ the
relation (1) converts into

$\displaystyle\mathcal{W}(X|\mu,\nu)=\frac{1}{2\pi|\nu|}\Big{|}\int\limits_{-\infty}^{\infty}\Psi(y)\exp{(\frac{i\mu}{2\nu}y^{2}-\frac{iX}{\nu}y)}dy\Big{|}^{2}.$
(6)

The representation (5) is connected with the Weyl and star-product
quantization [35]. It is known that the Weyl symbol of the density matrix is
explicitly the Wigner function. It is shown in [36] that the symplectic
tomogram is related to the quantum state expressed in terms of its Wigner
function $W(q,p)$ with an integral transform. The parameters $\mu$ and $\nu$
describe an ensemble of rotated and scaled reference frames, in which the
observable $X$ is measured. For $\mu=\cos{\phi}$ and $\nu=\sin{\phi}$, the pdf
$\mathcal{W}(X|\mu,\nu)$ is the distribution for homodyne output variable used
in optical tomography [37, 38]. The procedure of balanced homodyne photon
detection is based on mixing of a measurable (weak) field and a strong
coherent field with varying phase $\phi$ on the beam splitter. In this case,
the measurable observable is $\hat{X}=\hat{q}\cos{\phi}+\hat{p}\sin{\phi}$.
The angle $\phi$ could be interpreted as a rotation angle of the phase space.
Relations (1)-(6) for symplectic tomogram are transformed to equivalent
relations for optical tomograms. Note that the symplectic tomogram is a
function of two parameters $(\mu,\nu)$ of the $Sp(2,\mathrm{R})$ group
parameterization, and the optical tomogram is a function of the parameter
$\phi$.

Alternatively, the parameters $\mu$ and $\nu$ can be expressed in the form
$s\cos{\phi}$, $s^{-1}\sin{\phi}$, where $s>0$ is a real squeezing parameter
and $\phi$ is a rotation angle. Then the variable $X$ is identical to the
position measured in the new reference frame in the phase-space with axis $sq$
and $s^{-1}p$ and after the scaling the axis are rotated by an angle $\phi$.
Thus the tomogram implies the pdf of the random position $X$ measured in the
new (scaled and rotated) reference frame in the phase-space [39].

The process of inverting the raw observed data in order to arrive at a form of
the quantum state is both delicate and intriguing. Deterministic and
nondeterministic inversion methods are the two primary categories of inversion
techniques. In the deterministic methods an experimentally determined pdf
(tomogram) is used to determine the matrix elements of the density matrix, by
a direct mathematical inversion of (5). However this strongly depends on how
well we estimated the pdf based on the measurement of $X$ for varying $\phi$.
The nondeterministic technique aims to directly estimate the quantum state, as
opposed to using classical distributions as an intermediate object [40, 41].

##  3 Characteristic Function

Since the tomogram is the pdf of the variable $X$ depending on real parameters
$\mu$ and $\nu$, we can introduce its characteristic function. According to
the definition the characteristic function

$\displaystyle\phi_{X}(t)\equiv\int\limits_{-\infty}^{\infty}\mathcal{W}(X)e^{itX}dX$
(7)

is the Fourier transform of the pdf. It is known that a random variable
$X_{n}$ weakly converges to a random variable $X$ if and only if, for any $t$
the characteristic function $\phi_{X_{n}}(t)$ converges to the characteristic
function $\phi_{X}(t)$ [42]. Further, we omit $X$ in the notation
$\phi_{X}(t)$.

The characteristic function of any real-valued random variable completely
defines its probability distribution. According to definition, it is non-
vanishing in a region around zero ($\phi(0)=1$) and it is bounded
$|\phi(t)|\leq 1$, $\forall t$. We can introduce the characteristic function
$\phi(1;\mu,\nu)$ at point $t=1$ corresponding to the tomogram pdf (1):

$\displaystyle\phi(1;\mu,\nu)\equiv\int\limits_{-\infty}^{\infty}\mathcal{W}(X|\mu,\nu)e^{iX}dX.$
(8)

Then we can rewrite the density operator (5) as follows:

$\displaystyle\hat{\rho}=\frac{1}{2\pi}\int\phi(1;\mu,\nu)e^{-i(\mu\hat{q}+\nu\hat{p})}d\mu
d\nu,$ (9)

and, taking the trace from both sides, we can conclude

$\displaystyle\text{Tr}{\hat{\rho}}=\int
e^{i\frac{\mu\nu}{2}}\phi(1;\mu,\nu)\delta(\mu)\delta(\nu)d\mu
d\nu=\phi(1;0,0),$ (10)

where we used $Tr(e^{-i(\nu\hat{p}+\mu\hat{q})})=2\pi
e^{i\frac{\mu\nu}{2}}\delta(\mu)\delta(\nu)$. Thereby, the trace of the
density matrix of a quantum state is the characteristic function with
$\mu=0,\nu=0$. That gives us the condition on the characteristic function of
the quantum state: $\phi(1;0,0)=1$.

Using (5), we can write the product of two density matrices, $\hat{\rho}_{1}$
and $\hat{\rho}_{2}$ corresponding to two quantum states. Taking the trace, we
get

$\displaystyle\text{Tr}{(\hat{\rho}_{1}\hat{\rho}_{2})}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\iint\phi_{1}(1;\mu_{1},\nu_{1})\phi_{2}(1;-\mu_{1},-\nu_{1})d\mu_{1}d\nu_{1},$
(11)

where $\phi_{1}(1;\mu_{1},\nu_{1})$ and $\phi_{2}(1;\mu_{1},\nu_{1})$ are the
characteristic functions, corresponding to the tomograms of the states.

Since the modulus of the characteristic function is bounded by one, we can be
sure that the latter integral is also always bounded by one that fully
coincides with the left-hand side trace upper bound. However, the integral can
be negative since the characteristic function in general can take negative
values. The left-hand side trace from the product of two density matrices is
always non-negative, so the condition on the characteristic function of the
quantum state holds:

$\displaystyle
0\leq\frac{1}{2\pi}\iint\limits_{-\infty}^{\infty}\phi_{1}(1;\mu,\nu)\phi_{2}(1;-\mu,-\nu)d\mu
d\nu\leq 1.$ (12)

Formula (11) is an analogy of the density matrix distinguishability test.

######  Theorem 3.1.

Two quantum states with the density matrices $\hat{\rho}_{1}$ and
$\hat{\rho}_{2}$ are distinguishable if the integral from their characteristic
functions (11) is different from one.

When $\hat{\rho}_{1}=\hat{\rho}_{2}$, holds, the left-hand side of (11) is the
purity parameter

$\displaystyle\text{Tr}{(\hat{\rho}^{2})}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\iint\limits_{-\infty}^{\infty}\phi(1;\mu,\nu)\phi(1;-\mu,-\nu)d\mu
d\nu.$ (13)

Hence, if one knows the characteristic functions of the states, one can
calculate the purity or trace product of the density matrices.

The fidelity of two quantum states can be written in the same way as:

$\displaystyle\omega_{\psi_{1}\rightarrow\psi_{2}}=\text{Tr}{(\rho_{\phi_{2}}\rho_{\psi_{1}}^{\dagger})}=\frac{1}{2\pi}\int\phi_{2}(1;\mu,\nu)\phi_{1}(-1;-\mu,-\nu)d\mu
d\nu.$ (14)

Since the density matrix is hermitian, its matrix elements satisfy
$\rho(y,y^{\prime})=\rho^{\star}(y^{\prime},y)$. The matrix element of (5) is
(see the deduction in Appendix A)

$\displaystyle\rho(y,y^{\prime})=\frac{1}{2\pi}\int\phi(1;\mu,y-y^{\prime})\exp{\left(-i\frac{\mu(y+y^{\prime})}{2}\right)}d\mu.$
(15)

We can conclude that for a quantum state, the characteristic function must
satisfy:

$\displaystyle\phi(1;\mu,y-y^{\prime})=\phi(-1;\mu,y^{\prime}-y),\quad\forall
y,y^{\prime}.$ (16)

The diagonal elements of the density matrix can be written as

$\displaystyle\rho(y,y)=\frac{1}{2\pi}\int\phi(1;\mu,0)\exp{\left(-i\mu
y\right)}d\mu\equiv\hat{\phi}(1;y,0),$ (17)

where $\hat{\phi}(1;y,0)$ is a Fourier transform of $\phi(1;\mu,0)$. We can
conclude that $\hat{\phi}(1;y,0)\geq 0$.

We can summarise the latter results in the following theorem:

######  Theorem 3.2.

The characteristic function of a quantum state defines a density operator
$\hat{\rho}$ if it satisfies the following conditions:

$\displaystyle\text{Tr}{\hat{\rho}}=1:\quad\phi(1;0,0)=1,$ (18) $\displaystyle
0\leq\text{Tr}{\hat{\rho}^{2}}\leq 1:\quad
0\leq\quad\frac{1}{2\pi}\iint\limits_{-\infty}^{\infty}\phi(1;\mu,\nu)\phi(1;-\mu,-\nu)d\mu
d\nu\leq 1,$ (19)
$\displaystyle{\rho}^{\star}(y,y^{\prime})={\rho}(y^{\prime},y):\quad\phi(1;\mu,y-y^{\prime})=\phi(-1;\mu,y^{\prime}-y),\quad\forall
y,y^{\prime},$ (20) $\displaystyle\rho(y,y)\geq 0:\quad\hat{\phi}(1;y,0)\geq
0,\quad\forall y.$ (21)

For two density matrices $\hat{\rho}_{1}$, $\hat{\rho}_{2}$ the following
condition

$\displaystyle 0\leq\text{Tr}{(\hat{\rho}_{1}\hat{\rho}_{2})}\leq 1:\quad
0\leq\frac{1}{2\pi}\iint\limits_{-\infty}^{\infty}\phi_{1}(1;\mu,\nu)\phi_{2}(1;-\mu,-\nu)d\mu
d\nu\leq 1,$ (22)

holds.

Hence, if we want to reconstruct a state from a given tomographic pdf, we need
to ensure that its characteristic function satisfies the Theorem 3.2. Next, we
consider several examples.

###  3.1 Harmonic oscillator

The ground state of the harmonic oscillator (HO) is defined by the wave
function

$\displaystyle\Psi_{0}(y)=\pi^{-1/4}\exp{(-y^{2}/2)}.$ (23)

The corresponding tomogram can be easily calculated

$\displaystyle\mathcal{W}_{0}(X|\mu,\nu)=\frac{1}{\sqrt{\pi(\mu^{2}+\nu^{2})}}\exp{\left[-\frac{X^{2}}{\mu^{2}+\nu^{2}}\right]}.$
(24)

The latter is a Gaussian pdf $N(0,\frac{\mu^{2}+\nu^{2}}{2})$. The
corresponding characteristic function is

$\displaystyle{\phi}_{0}(1;\mu,\nu)=\exp{\left[-\frac{\mu^{2}+\nu^{2}}{4}\right]},$
(25)

and it is resistant to the change of the sign of $\mu$, $\nu$ of the
parameters. It is easy to check that the Theorem 3.2 is fulfilled.

The wave function of the excited state of the harmonic oscillator is

$\displaystyle\Psi_{n}(y)=\frac{1}{\sqrt[4]{\pi}\sqrt{2^{n}n!}}e^{-y^{2}/2}H_{n}(y),$
(26)

where $H_{n}(y)$ is the Hermite polynomial. We use the integral

$\displaystyle\int\limits_{-\infty}^{\infty}e^{-p(x-y)^{2}}H_{n}(cx)dx=\frac{\sqrt{\pi}(p-c^{2})^{n/2}}{p^{(n+1)/2}}H_{n}\left(cy\sqrt{\frac{p}{p-c^{2}}}\right)$
(27)

to find the tomogram of the excited oscillator state:

$\displaystyle\mathcal{W}_{n}(X|\mu,\nu)=\mathcal{W}_{0}(X|\mu,\nu)\frac{1}{2^{n}n!}\left(H_{n}\left(\frac{X}{\sqrt{\mu^{2}+\nu^{2}}}\right)\right)^{2}.$
(28)

To our best knowledge, the latter pdf is not arising in classical probability
theory in any context except quantum mechanics. In some rare sources, it is
named the Hermite-Gaussian pdf. We use the following integral:

$\displaystyle\int\limits_{-\infty}^{\infty}e^{-z^{2}+bz}H_{n}(z)H_{m}(z)dz=\sqrt{\pi}e^{\frac{b^{2}}{4}}n!m!\sum\limits_{k=0}^{\min{(n,m)}}\frac{2^{k}b^{n+m-2k}}{k!(n-k)!(m-k)!},$
(29)

to find its characteristic function

$\displaystyle{\phi}_{n}(1;\mu,\nu)$ $\displaystyle=$ $\displaystyle
e^{\frac{-(\mu^{2}+\nu^{2})}{4}}L_{n}\left(\frac{\mu^{2}+\nu^{2}}{2}\right),$
(30)

where we used the Laguerre polynomial series decomposition

$\displaystyle
L_{n}(x)=\sum\limits_{k=0}^{n}\frac{(-1)^{k}n!}{(k!)^{2}(n-k)!}x^{k}.$ (31)

One can check that the Theorem 3.2 is fulfilled.

##  4 Exponential family of pdfs

As we can see from the example of harmonic oscillators, the pdfs (24) and (28)
can be written in the following form:

$\displaystyle W(X|\eta)\equiv h(X)e^{\eta^{T}\tau(X)-A(\eta)},$ (32)

that is a general expression for an exponential family of pdfs for a given
vector of sufficient statistics $\tau(X)$ and normalization function $h(X)$.
For example, for the ground state of the harmonic oscillator, the functions
are the following: $\tau(X)=-X^{2}$, $\eta=-1/(\mu^{2}+\nu^{2})$,
$h(X)=1/\sqrt{\pi}$.

From the normalization condition, we can find the cumulant generating
function:

$\displaystyle e^{A(\eta)}=\int h(X)e^{\eta^{T}\tau(X)}dX.$ (33)

This shows that $A(\eta)$ is not a degree of freedom in the specification of
an exponential family density. It is determined once $\eta$, $\tau(X)$ and
$h(X)$ are determined:

$\displaystyle A(\eta)=\log{\left(\int h(X)e^{\eta^{T}\tau(X)}dX\right)}.$
(34)

The characteristic function (7) in the point $t=1$ for such family of pdfs can
be written as follows:

$\displaystyle\phi(1;\eta(\mu,\nu))=\int
e^{iX}W(X|\eta)dX=e^{-A(\eta(\mu,\nu))}\int
e^{iX+\eta(\mu,\nu)^{T}\tau(X)}h(X)dX,$ (35)

where $\eta(\mu,\nu)$ is a vector, depending on $\mu$, $\nu$. If one wants to
check the known pdf from the exponential family to generate the quantum state,
one needs to check that the latter characteristic function satisfies the
Theorem 3.2. We deduce the following:

######  Corollary 4.1.

The characteristic function from an exponential family of pdfs (35) describes
a quantum state defined by a density matrix $\hat{\rho}$ if:

$\displaystyle\int e^{iX+\eta(0,0)^{T}\tau(X)}h(X)dX=e^{A(\eta(0,0))},$ (36)
$\displaystyle 0$ $\displaystyle\leq$
$\displaystyle\frac{1}{2\pi}\iint\limits_{-\infty}^{\infty}e^{-A(\eta(\mu,\nu))-A(\eta(-\mu,-\nu))}\int
e^{iX_{1}+\eta(\mu,\nu)^{T}\tau(X_{1})}h(X_{1})dX_{1}$ $\displaystyle\times$
$\displaystyle\int e^{iX_{2}+\eta(-\mu,-\nu)^{T}\tau(X_{2})}h(X_{2})dX_{2}d\mu
d\nu\leq 1,$ $\displaystyle e^{A(\eta(\mu,y^{\prime}-y))}\int
e^{iX+\eta(\mu,y-y^{\prime})^{T}\tau(X)}h(X)dX$ $\displaystyle=$
$\displaystyle e^{A(\eta(\mu,y-y^{\prime}))}\int
e^{-iX+\eta(\mu,y^{\prime}-y)^{T}\tau(X)}h(X)dX,\quad\forall y,y^{\prime},$
$\displaystyle\frac{1}{2\pi}\int e^{-A(\eta(\mu,0))}\int
e^{iX+\eta(\mu,0)^{T}\tau(X)}h(X)e^{-i\mu y}dXd\mu\geq 0,\quad\forall y.$ (39)

Hence, any pdf from the exponential family can be checked on the latter
conditions. Further, we provide more detailed analyses for the most important
pdfs known in probability theory.

###  4.1 Special cases

Let us observe $\tau(X)=X$, $X\in\mathbb{R}^{+}$, $A(\eta)$ is the $\log{}$ of
Laplace transform of $h(X)$. The Laplace transform of $h(X)$ is defined as
follows:

$\displaystyle H(s)=\int\limits_{0}^{\infty}h(X)e^{-sX}dX,$ (40)

where $s$ is a complex frequency domain parameter $s=\sigma+i\omega$ with
$\sigma,\omega\in\mathbb{R}$. A necessary condition for existence of the
integral is that $h(X)$ must be locally integrable on $(0,\infty)$ and
$\eta<0$. Then the characteristic function in $t=1$ is

$\displaystyle\phi(1;\eta)=\left(\int\limits_{0}^{\infty}h(X_{1})e^{\eta
X_{1}}dX_{1}\right)^{-1}\left(\int\limits_{0}^{\infty}h(X_{2})e^{(i+\eta)X_{2}}dX_{2}\right).$
(41)

######  Example 4.1.

Let us observe $h(X)=X^{\alpha-1}$, $\alpha>0$ and $\eta=-p(\mu,\nu)$, where
$p(\mu,\nu)$ is a nonnegative function. Then the pdf is

$\displaystyle f(X|\mu,\nu)=X^{\alpha-1}e^{-p(\mu,\nu)X-A(p(\mu,\nu))},\quad
X\in\mathbb{R}^{+}.$ (42)

Then the characteristic function is

$\displaystyle\phi(1;\mu,\nu)=\frac{p^{\alpha}(\mu,\nu)}{(p(\mu,\nu)-i)^{\alpha}}.$
(43)

One can easily check that the latter type of functions does not satisfy the
first condition of the Theorem 3.2. Hence, this class of pdfs is not suitable
for a tomogram. This class includes important pdfs like exponential, gamma,
and $\chi^{2}$.

Let us observe $\tau(X)=(X,X^{2})^{T}$, $\eta=(\eta_{1},\eta_{2})$, $X\in R$.
Then the characteristic function is

$\displaystyle\phi(1;\eta)=\left(\int\limits_{-\infty}^{\infty}h(X_{1})e^{\eta_{1}X_{1}+\eta_{2}X_{1}^{2}}dX_{1}\right)^{-1}\left(\int\limits_{-\infty}^{\infty}h(X_{2})e^{(i+\eta_{1})X_{2}+\eta_{2}X^{2}_{2}}dX_{2}\right).$
(44)

######  Example 4.2.

For a special case of $h(X)=C$, where $C$ is a constant and
$\eta_{1}=p_{1}(\mu,\nu)$, $\eta_{2}=-p_{2}(\mu,\nu)$, where $p_{i}(\mu,\nu)$,
$i=1,2$ are nonnegative functions, the characteristic function is the
following:

$\displaystyle\phi(1;p(\mu,\nu))=e^{-\frac{2ip_{1}(\mu,\nu)-1}{4p_{2}(\mu,\nu)}}.$
(45)

If we take $p_{1}(\mu,\nu)=0$, $p_{2}(\mu,\nu)=-1/(\mu^{2}+\nu^{2})$ the case
corresponds to harmonic oscillator (25) and the tomogram is the Gaussian
distribution.

##  5 Pseudoharmonic Oscillator

The straight motion along the positive $x-$ semi-axis can be described by
potential

$\displaystyle
V^{1D}(a,x)=D_{\omega}\left(\frac{x}{x_{\omega}}-\sqrt{a}\frac{x_{\omega}}{x}\right)^{2},\quad
0\leq x\leq\infty,$ (46)

where the dimensionless coefficient $a\geq 0$,
$D_{\omega}=\frac{\hbar\omega}{2}$, $\omega$ is the confining frequency and
$x_{\omega}=\frac{\hbar}{m\omega}$. At $a=0$, equation (46) describes the
geometry of a particle confined to the right-hand half of the harmonic
oscillator of the frequency $\omega$. The problem (46) is cited in literature
as pseudoharmonic oscillator (PHO) [34]. The 1D Schrödinger equation is

$\displaystyle\frac{\hbar^{2}}{2m}\frac{d^{2}}{dx^{2}}\Psi_{n}(a;x)+V^{1D}(a,x)\Psi_{n}(a;x)=E_{n}(a)\Psi_{n}(a;x).$
(47)

For the $a$-dependent potential (46) the energies and the wave function are

$\displaystyle E_{n}(a)=\hbar\omega(2n+1+b-\sqrt{a}),$ (48)
$\displaystyle\Psi_{n}(a;x)=\frac{1}{\sqrt{x_{\omega}}}\left[\frac{2n!}{\Gamma(n+b+1)}\right]^{\frac{1}{2}}\left(\frac{x}{x_{\omega}}\right)^{b+\frac{1}{2}}\exp{\left(-\frac{1}{2}\frac{x^{2}}{x^{2}_{\omega}}\right)}L_{n}^{(\eta)}\left(\frac{x^{2}}{x^{2}_{\omega}}\right),$

where $b=\frac{1}{2}\sqrt{1+4a}$, $\Gamma(z)$ is a $\Gamma$-function,
$L_{n}^{(b)}$ is an $n$th order associated Laguerre polynomial. The divergence
of the potential at the left edge $x=0$ forces the function to vanish there at
any $a$ and $n$: $\Psi_{n}(a;0)=0$. At $a=0$ the spatial dependence reads

$\displaystyle\Psi_{n}(0;x)=\frac{(-1)^{n}}{x_{\omega}^{\frac{1}{2}}}\frac{1}{2^{2n+\frac{1}{2}}}\left[\frac{1}{n!\Gamma(n+\frac{3}{2})}\right]^{\frac{1}{2}}\exp{\left(-\frac{1}{2}\frac{x^{2}}{x^{2}_{\omega}}\right)}H_{2n+1}\left(\frac{x}{x_{\omega}}\right).$
(49)

Then, one can find the tomogram $\mathcal{W}_{n}(X|a,\mu,\nu)$ of (49):

$\displaystyle\mathcal{W}_{n}(X|0,\mu,\nu)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\frac{1}{2\pi|\nu|}\frac{1}{x_{\omega}}\frac{1}{2^{4n+1}}\frac{1}{n!\Gamma(n+\frac{3}{2})}I(X),$
(50) $\displaystyle I(X)$ $\displaystyle\equiv$
$\displaystyle\Bigg{|}\int\limits_{0}^{\infty}\exp{\left[-x^{2}\left(\frac{1}{2x^{2}_{\omega}}-\frac{i\mu}{2\nu}\right)-\frac{iX}{\nu}x\right]}H_{2n+1}\left(\frac{x}{x_{\omega}}\right)dx\Bigg{|}^{2}.$

The closed form of the latter tomogram is derived in Appendix B:

Figure 1: The tomogram $\mathcal{W}_{n}(X|0,\mu,\nu)$ for the PHO (49) is
plotted for $n\in\\{0,1,2,10\\}$. The ground state tomogram is a Gaussian pdf
with the shifted center, while the excited states are the multi-mode pdfs from
the exponential family.

$\displaystyle\mathcal{W}_{n}(X|0,\mu,\nu)\\!\\!=\\!\\!\frac{1}{\pi}\frac{((2n+1)!)^{2}}{n!(n+\frac{1}{2})!}\frac{\nu^{2n+1}x_{\omega}}{(\nu^{2}+\mu^{2}x^{4}_{\omega})^{n+1}}e^{-\frac{X^{2}x^{2}_{\omega}}{2\left(\nu^{2}+\mu^{2}x^{4}_{\omega}\right)}}$
$\displaystyle\times$
$\displaystyle\Bigg{[}\sum_{m_{1},m_{2}=0}^{n}\frac{(-1)^{m_{1}+m_{2}}}{{m_{1}}!{m_{2}}!}\left(\frac{x^{2}_{\omega}}{16\nu}\right)^{{m_{1}}+m_{2}}\left(\nu-i\mu
x^{2}_{\omega}\right)^{m_{1}}\left(\nu+i\mu x^{2}_{\omega}\right)^{m_{2}}$
$\displaystyle\times$ $\displaystyle
D_{-(2n-2{m_{1}}+2)}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)D_{-(2n-2{m_{2}}+2)}\left(\frac{iXx_{\omega}\sqrt{\nu-i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)\Bigg{]}.$

This pdf is a, so called, normal-exponential-gamma distribution type. It is
depicted in Fig. 1 for the excitations $n=\\{0,1,2,10\\}$. The distribution is
not symmetrical with the center of the pdf shifted to the right side of the
$X$-axis, whereas the pdf of a standard harmonic oscillator is centered at
$X=0$. The characteristic function of the latter distribution defined by
transformation (8) is quite cumbersome. We have numerically verified that it
satisfies the conditions of Theorem 3.2.

Next we move to the tomogram for the general case (5). According to the
definition it can be written as

$\displaystyle\mathcal{W}(X|a,\mu,\nu)=\frac{1}{2\pi|\nu|}\frac{1}{x_{\omega}}\left[\frac{2n!}{\Gamma(n+b+1)}\right]\Bigg{|}\int\limits_{0}^{\infty}\left(\frac{y}{x_{\omega}}\right)^{b+\frac{1}{2}}e^{\left(-\frac{1}{2}\frac{y^{2}}{x^{2}_{\omega}}\right)}L_{n}^{(b)}\left(\frac{y^{2}}{x^{2}_{\omega}}\right)e^{(\frac{i\mu}{2\nu}y^{2}-\frac{iX}{\nu}y)}dy\Bigg{|}^{2}.$
(52)

a)

b)

Figure 2: Comparison of the tomograms of HO and the PHO for
$a\in\\{0,10,100,1000\\}$, a) The ground state ($n=0$) tomograms are the
Gaussian pdfs, where for the PHO the center of the pdf is shifted. b) the
first excited state ($n=1$) tomograms are from the exponential family of pdf
with multiple modes.

The closed form is

$\displaystyle\mathcal{W}_{n}(X|a,\mu,\nu)=\frac{x_{\omega}}{2\pi|\nu|}\left[\frac{2n!}{\Gamma(n+b+1)}\right]\exp{\left(-\frac{X^{2}x^{2}_{\omega}}{2(\nu^{2}+\mu^{2}x^{4}_{\omega})}\right)}\left(\frac{(b+1)_{n}}{n}\right)^{2}$
$\displaystyle\times$
$\displaystyle\Bigg{|}\sum\limits_{k=0}^{\infty}\frac{(n)_{k}}{(b+1)_{k}k!}\Gamma\left(2k+b+\frac{3}{2}\right)\left(\frac{\nu-i\mu
x^{2}_{\omega}}{\nu}\right)^{-(k+\frac{b}{2}+\frac{3}{4})}D_{-(2k+b+\frac{3}{2})}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)\Bigg{|}^{2}.$

Using the latter representation, we can examine the dependence of
$\mathcal{W}_{n}(X|a,\mu,\nu)$ from $a$ and compare it to the
$\mathcal{W}_{n}(X|\mu,\nu)$ pdf. Both functions are depicted in Fig. 2 for
$n=\\{0,1\\}$. We can observe that the larger the value of $a$ becomes, the
further the function $\mathcal{W}_{n}(X|a,\mu,\nu)$ is shifted to the right.
The shape of the pdf changes slightly with the growth of $a$. However, for the
first exited state, the mass of the picks of the distribution is slightly
changing.

We observe that, both for HO and PHO, the tomograms deviate from Gaussian
behavior even for the first excited state. The distribution becomes multi-
modal and challenging to estimate. Subsequently, we will examine the case of
an even more complex distribution arising in superposition states.

##  6 Crystallized cat states

The even and odd coherent states are introduced in [43] from the Glauber
coherent states $\ket{\alpha}$ of the harmonic oscillator as:

$\displaystyle\ket{\alpha_{\pm}}=N_{\pm}(\ket{\alpha}\pm\ket{-\alpha}),\quad
N_{\pm}=\left(2(1\pm e^{-2|\alpha|^{2}})\right)^{-1/2}.$ (54)

Their generalization is provided in [44], using the Abelian symmetry group
$C_{3}$ with three rotation elements $\\{1,e^{2\pi i/3},e^{4\pi i/3}\\}$,
acting on the coherent states, giving the following state:

$\displaystyle\ket{\psi}=N\sum\limits_{j=1}^{3}\ket{\psi_{j}},\quad\ket{\psi_{1}}=\ket{\alpha},\quad\ket{\psi_{2}}=\ket{\alpha
e^{2\pi i/3}},\quad\ket{\psi_{3}}=\ket{\alpha e^{4\pi i/3}}.$ (55)

The coordinate representation of the coherent state reads as

$\displaystyle\psi_{\alpha}(x)=\langle
x|\alpha\rangle=\pi^{-1/4}\exp{\left(-\frac{x^{2}}{2}-\frac{|\alpha|^{2}}{2}+\sqrt{2}\alpha
x-\frac{\alpha^{2}}{2}\right)},$ (56)

and its tomogram is known to be equal to [26]

$\displaystyle\mathcal{W}(X|{\alpha},\mu,\nu)=\frac{e^{-|\alpha|^{2}}e^{\frac{(\nu+i\mu)^{2}\alpha^{2}+(\nu-i\mu)^{2}(\alpha^{\star})^{2}}{2(\nu^{2}+\mu^{2})}}}{\sqrt{\pi(\nu^{2}+\mu^{2})}}e^{-\frac{X^{2}}{\nu^{2}+\mu^{2}}}e^{\frac{\sqrt{2}iX\left((\nu-i\mu)\alpha^{\star}-(\nu+i\mu)\alpha\right)}{\nu^{2}+\mu^{2}}}.$
(57)

The characteristic function is

$\displaystyle{\phi}(1;{\alpha},\mu,\nu)=e^{-\frac{\nu^{2}+\mu^{2}}{4}}e^{-\frac{(\nu-i\mu)\alpha^{\star}-(\nu+i\mu)\alpha}{\sqrt{2}}}.$
(58)

We can write every state in (55) as

$\displaystyle\psi_{j}(x)=N_{j}\exp{\left(-A_{j}x^{2}+B_{j}x+C_{j}\right)},\quad
j=1,2,3,$ (59) $\displaystyle A_{j}=\frac{1}{2},\quad B_{j}=\sqrt{2}\alpha
e^{\frac{2\pi i}{3}(j-1)},\quad
C_{j}=-\frac{|\alpha|^{2}}{2}-\frac{\alpha^{2}e^{\frac{4\pi
i}{3}(j-1)}}{2},\quad N_{j}=\pi^{-1/4}.$

According to (6) the tomogram can be written as

$\displaystyle\mathcal{W}_{ccat}(X|\alpha,\mu,\nu)$ $\displaystyle=$
$\displaystyle\frac{|N|^{2}e^{-|\alpha|^{2}}}{\sqrt{\pi(\nu^{2}+\mu^{2})}}e^{-\frac{X^{2}}{\nu^{2}+\mu^{2}}}$
$\displaystyle\times$
$\displaystyle\sum\limits_{j,k=1}^{3}\exp{\left(\frac{\sqrt{2}iX\left((\nu-i\mu)\alpha^{\star}e^{-\frac{2\pi
i}{3}(k-1)}-(\nu+i\mu)\alpha e^{\frac{2\pi
i}{3}(j-1)}\right)}{\nu^{2}+\mu^{2}}\right)}$ $\displaystyle\times$
$\displaystyle\exp{\left(\frac{\left((\nu+i\mu)^{2}\alpha^{2}e^{\frac{4\pi
i}{3}(j-1)}+(\nu-i\mu)^{2}(\alpha^{\star})^{2}e^{-\frac{4\pi
i}{3}(k-1)}\right)}{2(\nu^{2}+\mu^{2})}\right)}.$

One can see that if $k,j=1$ the tomogram coincides with the single coherent
state tomogram (57). The pdf (6) is a mixture of the Gaussian distributions
[45]. Note that the sum of Gaussian pdfs in general is not a Gaussian pdf
itself.

The characteristic function is the following:

$\displaystyle\phi_{ccat}(1;\alpha,\mu,\nu)$ $\displaystyle=$
$\displaystyle|N|^{2}e^{-\frac{\nu^{2}+\mu^{2}}{4}}\sum\limits_{j,k=1}^{3}e^{-|\alpha|^{2}(\exp{(-\frac{2\pi
i}{3}(k-j)}+1)}$ $\displaystyle\times$ $\displaystyle
e^{-\frac{(\nu-i\mu)\alpha^{\star}\exp{(-\frac{2\pi
i}{3}(k-1))}-(\nu+i\mu)\alpha\exp{(\frac{2\pi i}{3}(j-1))}}{\sqrt{2}}}.$

One can check that it satisfies the Theorem 3.2.

The Gottesman-Kitaev-Preskill (GKP) code is specifically designed for
continuous-variable quantum systems, and its states are characterized by
coherent superpositions [46]. In particular, the GKP code employs states that
are approximate eigenstates of position and momentum operators. The GKP states
are typically defined as Gaussian wavepacket superpositions, and they serve as
an encoding for quantum information. Since the GKP code space is an abstract
construction and the code words are not normalised, in practise the coherent
state approximation is used to prepare the state lying entirely in the GKP
code space [47]. Thus, the tomogram of a specific superposition of coherent
states serves as an approximation to the tomogram of the theoretical GKP
state. We see in this a potential extension of the article into the realm of
error correction codes.

| Quantum
---
state
| Tomogram (pdf)
---
Characteristic function
HO (ground state) |  $\mathcal{W}_{0}(X|\mu,\nu)=\frac{1}{\sqrt{\pi(\mu^{2}+\nu^{2})}}\exp{\left[-\frac{X^{2}}{\mu^{2}+\nu^{2}}\right]}$,
| ${\phi}_{0}(1|\mu,\nu)=\exp{\left[-\frac{\mu^{2}+\nu^{2}}{4}\right]}$
HO (excited state) | $\mathcal{W}_{n}(X|\mu,\nu)=\mathcal{W}_{0}(X|\mu,\nu)\frac{1}{2^{n}n!}\Big{(}H_{n}\left(\frac{X}{\sqrt{\mu^{2}+\nu^{2}}}\right)\Big{)}^{2}$
|
${\phi}_{n}(1;\mu,\nu)={\phi}_{0}(1|\mu,\nu)L_{n}\left(\frac{\mu^{2}+\nu^{2}}{2}\right)$
PHO ($a=0$) | $\mathcal{W}_{n}(X|\frac{1}{2},\mu,\nu)\\!\\!=\\!\\!\frac{1}{\pi}\frac{((2n+1)!)^{2}}{n!(n+\frac{1}{2})!}\frac{\nu^{2n+1}x_{\omega}}{(\nu^{2}+\mu^{2}x^{4}_{\omega})^{n+1}}e^{-\frac{X^{2}x^{2}_{\omega}}{2\left(\nu^{2}+\mu^{2}x^{4}_{\omega}\right)}}$
|
$\times\Bigg{[}\sum_{m_{1},m_{2}=0}^{n}\frac{(-1)^{m_{1}+m_{2}}}{{m_{1}}!{m_{2}}!}\left(\frac{x^{2}_{\omega}}{16\nu}\right)^{{m_{1}}+m_{2}}\left(\nu-i\mu
x^{2}_{\omega}\right)^{m_{1}}\left(\nu+i\mu x^{2}_{\omega}\right)^{m_{2}}$
| $\times D_{-(2n-2{m_{1}}+2)}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)D_{-(2n-2{m_{2}}+2)}\left(\frac{iXx_{\omega}\sqrt{\nu-i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)\Bigg{]}.$
PHO ($\forall a$) | $\mathcal{W}_{n}(X|\eta,\mu,\nu)=\frac{x_{\omega}}{2\pi|\nu|}\left[\frac{2n!}{\Gamma(n+\eta+1)}\right]\exp{\left(-\frac{X^{2}x^{2}_{\omega}}{2(\nu^{2}+\mu^{2}x^{4}_{\omega})}\right)}\left(\frac{(\eta+1)_{n}}{n}\right)^{2}$
|
$\Bigg{|}\sum\limits_{k=0}^{\infty}\frac{(n)_{k}}{(\eta+1)_{k}k!}\Gamma\left(2k+\eta+\frac{3}{2}\right)\left(\frac{\nu-i\mu
x^{2}_{\omega}}{\nu}\right)^{-(k+\frac{\eta}{2}+\frac{3}{4})}\\!\\!\\!D_{-(2k+\eta+\frac{3}{2})}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)\Bigg{|}^{2}$
Coherent state | $\mathcal{W}(X|{\alpha},\mu,\nu)=\mathcal{W}_{0}(X|\mu,\nu)e^{-|\alpha|^{2}}e^{\frac{(\nu+i\mu)^{2}\alpha^{2}+(\nu-i\mu)^{2}(\alpha^{\star})^{2}}{2(\nu^{2}+\mu^{2})}}e^{\frac{\sqrt{2}iX\left((\nu-i\mu)\alpha^{\star}-(\nu+i\mu)\alpha\right)}{\nu^{2}+\mu^{2}}}$
|
${\phi}(1;{\alpha},\mu,\nu)={\phi}_{0}(1;\mu,\nu)e^{-\frac{(\nu-i\mu)\alpha^{\star}-(\nu+i\mu)\alpha}{\sqrt{2}}}$
Crystallized | $\mathcal{W}_{ccat}(X|\alpha,\mu,\nu)=\mathcal{W}_{0}(X|\mu,\nu)e^{-|\alpha|^{2}}|N|^{2}$
cat states | $\\!\\!\times\\!\\!\sum\limits_{j,k=1}^{3}e^{\frac{(\nu+i\mu)^{2}\alpha^{2}e^{\frac{4\pi i}{3}(j-1)}+(\nu-i\mu)^{2}(\alpha^{\star})^{2}e^{-\frac{4\pi i}{3}(k-1)}}{2(\nu^{2}+\mu^{2})}}\\!\\!e^{\frac{\sqrt{2}iX((\nu-i\mu)\alpha^{\star}e^{-\frac{2\pi i}{3}(k-1)}-(\nu+i\mu)\alpha e^{\frac{2\pi i}{3}(j-1)})}{\nu^{2}+\mu^{2}}}$
| ${\phi}_{ccat}(1;{\alpha},\mu,\nu)=|N|^{2}{\phi}_{0}(1;\mu,\nu)$
| $\sum\limits_{j,k=1}^{3}e^{-|\alpha|^{2}(\exp{(-\frac{2\pi
i}{3}(k-j)}+1)}e^{-\frac{(\nu-i\mu)\alpha^{\star}\exp{(-\frac{2\pi
i}{3}(k-1))}-(\nu+i\mu)\alpha\exp{(\frac{2\pi i}{3}(j-1))}}{\sqrt{2}}}$
Table 1: Tomograms and characteristic functions for the ground and excited
states of the oscillator and pseudoharmonic oscillator, as well as for
coherent states and superpositions of cat states.

##  7 Summary

For the accurate reconstruction of a quantum state, various techniques are
employed. The tomographic pdf, being the true distribution function,
completely characterizes the state and have potential to be a convenient
method for quantum state reconstruction. When applying classical probability
theory methods for estimating the pdf, it is often necessary to know which
distribution family we are working with. In this article, we thoroughly
investigated the conditions (Theorem 3.2) that the characteristic function of
the distribution must satisfy to describe a physical system. This
characteristic function is obtained as the Fourier transform of the pdf
tomogram, aiming to describe a physical quantum state. These conditions
strictly limit the family of pdfs suitable for describing quantum states. We
specifically explored the exponential family of densities, introducing
conditions that the characteristic function of a quantum state from this
family must satisfy. We tested several well-known distributions, such as
exponential, gamma, and $\chi^{2}$, demonstrating their unsuitability for
generating quantum states.

Using the characteristic function, it is possible to express key features of a
quantum state, e.g., the purity parameter. Additionally, we introduced a trace
distance between two states based on their characteristic functions.
Consequently, knowledge of the estimate of the characteristic function can
serve as a test for the purity or proximity of two quantum states.

We investigated well-known states such as the quantum oscillator (symmetric
and asymmetric), coherent states, cat states, and their superpositions. The
tomograms of these states are associated with Gaussian and Hermite-Gaussian
distributions or mixtures of these distributions, as shown in Table 1. Their
characteristic functions satisfy the Theorem 3.2 as expected.

The study clearly demonstrates that the family of densities capable of
generating quantum states is extremely limited. Further development of the
article could involve a detailed examination of well-known density
distributions, exploring their potential as tomograms for quantum states. Our
intuition suggests that the presence of a quadratic exponential term (as in
the tomography of the oscillator and all other tomograms in Table 1) is likely
mandatory. It would be particularly interesting to investigate the power-law
family of distributions, which finds extensive application in classical
mathematics too. This study provides insights into the family of distributions
encountered in quantum mechanics. With this information, one can employ both
parametric and non-parametric methods to estimate the tomograms of quantum
states. This investigation will be the focus of our future work.

##  8 Acknowledgments

L.M. was partly supported by the Netherlands Organisation for Scientific
Research (NWO/OCW), as part of the Quantum Software Consortium program
(project number 024.003.037 / 3368). This research work was partly supported
by the Roadmap for the Development of Quantum Technologies in Russian
Federation, contract No. 868-1.3-15/15-2021.

## References

  * [1] Erwin Schrödinger.  Quantisierung als eigenwertproblem.  Uspekhi Fizicheskikh Nauk, 122(4):621–632, 1977.
  * [2] L. Landau.  Das dämpfungsproblem in der wellenmechanik.  Zeitschrift für Physik, 45(5-6):430–441, 1927.
  * [3] J. Neumann.  Wahrscheinlichkeitstheoretischer aufbau der quantenmechanik.  Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, 1927:245–272, 1927.
  * [4] P. Dirac.  The principles of quantum mechanics.  Number 27. Oxford university press, 1981.
  * [5] B.W. Silverman.  Density estimation for statistics and data analysis, volume 26.  CRC press, 1986.
  * [6] E. Wigner.  On the quantum correction for thermodynamic equilibrium.  Phys. Rev., 40:749–759, Jun 1932.
  * [7] K. Husimi.  Some formal properties of the density matrix.  Proceedings of the Physico-Mathematical Society of Japan. 3rd Series, 22(4):264–314, 1940.
  * [8] K. Yutaka.  A new phase-space distribution function in the statistical theory of the electromagnetic field.  J. Math. Phys., 6(12):1913–1915, 1965.
  * [9] T. Linowski and Ł. Rudnicki.  Relating the glauber-sudarshan, wigner and husimi quasiprobability distributions operationally through the quantum limited amplifier and attenuator channels.  arXiv preprint arXiv:2312.08216, 2023.
  * [10] R.J. Glauber.  Photon correlations.  Phys. Rev. Lett., 10:84–86, Feb 1963.
  * [11] E.C.G. Sudarshan.  Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams.  Phys. Rev. Lett., 10:277–279, Apr 1963.
  * [12] J.G. Kirkwood.  Quantum statistics of almost classical assemblies.  Phys. Rev., 44:31–37, Jul 1933.
  * [13] L. Praxmeyer and K. Wodkiewicz.  Quantum interference in the Kirkwood-Rihaczek representation.  Opt. Comm., 223(4):349–365, 2003.
  * [14] H. Margenau and R.N. Hill.  Correlation between measurements in quantum theory.  Prog. Theor. Phys., 26(5):722–738, 11 1961.
  * [15] C.H. Page.  Instantaneous power spectra.  J. Appl. Phys, 23(1):103–106, 1952.
  * [16] H. Choi and W.J. Williams.  Improved time-frequency representation of multicomponent signals using exponential kernels.  IEEE Trans. Acoust. Speech Signal Process., 37:862–871, 1989.
  * [17] L. Cohen.  Time-frequency distributions-a review.  Proceedings of the IEEE, 77(7):941–981, 1989.
  * [18] L. Cohen.  Generalized phase-space distribution functions.  J. Math. Phys., 7(5):781–786, 1966.
  * [19] S. Mancini, V.I. Man’ko, and P. Tombesi.  Symplectic tomography as classical approach to quantum systems.  Phys. Lett. A, 213(1):1 – 6, 1996.
  * [20] S. Mancini, V.I. Man’ko, and P. Tombesi.  Wigner function and probability distribution for shifted and squeezed quadratures.  Quantum and Semiclassical Optics: Journal of the European Optical Society Part B, 7(4):615–623, aug 1995.
  * [21] V.I. Man’ko and L.A. Markovich.  Integral transforms between tomogram and quasi-probability functions based on quantizer-dequantizer operators formalism.  Journal of Mathematical Physics, 61(10), 2020.
  * [22] A Ibort, V I Man’ko, G Marmo, A Simoni, and F Ventriglia.  An introduction to the tomographic picture of quantum mechanics.  Physica Scripta, 79(6):065013, jun 2009.
  * [23] J. Řeháček, D. Mogilevtsev, and Z. Hradil.  Tomography for quantum diagnostics.  New Journal of Physics, 10(4):043022, apr 2008.
  * [24] N. Fabre, A. B Klimov, G. Leuchs, and L.L. Sánchez-Soto.  Local sampling of the su (1, 1) wigner function.  AVS Quantum Science, 5(1), 2023.
  * [25] Mattia Walschaers.  Non-gaussian quantum states and where to find them.  PRX Quantum, 2:030204, Sep 2021.
  * [26] O. Man’ko and V.I. Man’ko.  Quantum states in probability representation and tomography.  Journal of Russian Laser Research, 18(5):407–444, 1997.
  * [27] O. Barndorff-Nielsen.  Information and exponential families: in statistical theory.  John Wiley & Sons, 2014.
  * [28] O.V. Man’ko and V.I. Man’ko.  Inverted oscillator quantum states in the probability representation.  Entropy, 25(2):217, 2023.
  * [29] M.A. Man’ko and V.I. Man’ko.  Quantum oscillator at temperature t and the evolution of a charged-particle state in the electric field in the probability representation of quantum mechanics.  Entropy, 25(2):213, 2023.
  * [30] Jesus R. Moya-Cessa, Luis Raul Berriel-Valdos, and Héctor Manuel Moya-Cessa.  Optical production of the husimi function of two gaussian functions.  Applied Mathematics and Information Sciences, 2:309–316, 2008.
  * [31] B. Yurke and D. Stoler.  Generating quantum mechanical superpositions of macroscopically distinguishable states via amplitude dispersion.  Phys. Rev. Lett., 57:13–16, Jul 1986.
  * [32] Daniel Gottesman, Alexei Kitaev, and John Preskill.  Encoding a qubit in an oscillator.  Phys. Rev. A, 64:012310, Jun 2001.
  * [33] J.A. López-Saldívar, M.A. Man’ko, and V.I. Man’ko.  Symplectic tomographic probability distribution of crystallized schrödinger cat states.  Physics Letters A, 434:128044, 2022.
  * [34] O. Olendski.  One-dimensional pseudoharmonic oscillator: classical remarks and quantum-information theory.  Journal of Physics Communications, 7(4):045002, 2023.
  * [35] A. Ibort, V.I. Man’Ko, G. Marmo, A. Simoni, and F. Ventriglia.  An introduction to the tomographic picture of quantum mechanics.  Physica Scripta, 79(6):065013, 2009.
  * [36] S. Mancini, V.I. Man’ko, and P. Tombesi.  Wigner function and probability distribution for shifted and squeezed quadratures.  Quantum and Semiclassical Optics: Journal of the European Optical Society Part B, 7(4):615, 1995.
  * [37] K. Vogel and H. Risken.  Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase.  Phys. Rev. A, 40:2847–2849, Sep 1989.
  * [38] U. Leonhardt and H. Paul.  High-accuracy optical homodyne detection with low-efficiency detectors:” preamplification” from antisqueezing.  Physical review letters, 72(26):4086, 1994.
  * [39] A. Del Campo, V.I. Man’ko, and G. Marmo.  Symplectic tomography of ultracold gases in tight waveguides.  Physical Review A, 78(2):025602, 2008.
  * [40] K.R.W. Jones.  Fundamental limits upon the measurement of state vectors.  Physical Review A, 50(5):3682, 1994.
  * [41] R. Schack, T.A. Brun, and C.M. Caves.  Quantum bayes rule.  Physical Review A, 64(1):014305, 2001.
  * [42] A. N. Kolmogorov and A.T. Bharucha-Reid.  Foundations of the theory of probability: Second English Edition.  Courier Dover Publications, 2018.
  * [43] V.V. Dodonov, I.A. Malkin, and V.I. Man’Ko.  Even and odd coherent states and excitations of a singular oscillator.  Physica, 72(3):597–615, 1974.
  * [44] V.V. Dodonov, I.A. Malkin, and V.I. Man’ko.  Even and odd coherent states and excitations of a singular oscillator.  Physica, 72(3):597–615, 1974.
  * [45] CA Robertson and JG Fryer.  Some descriptive properties of normal mixtures.  Scandinavian Actuarial Journal, 1969(3-4):137–146, 1969.
  * [46] A.L. Grimsmo and S. Puri.  Quantum error correction with the gottesman-kitaev-preskill code.  PRX Quantum, 2(2):020101, 2021.
  * [47] N.C. Menicucci.  Fault-tolerant measurement-based quantum computing with continuous-variable cluster states.  Phys. Rev. Lett., 112:120504, Mar 2014.

##  Appendix A Matrix elements of the density matrix

To write the condition $\rho{(y,y^{\prime})}=\rho(y^{\prime},y)^{\star}$ in
terms of the characteristic function, we use the Baker- Campbell- Hausdorff
formula
$e^{ia\hat{q}}e^{ib\hat{p}}e^{\frac{iab\hbar}{2}}=e^{i\left(a\hat{q}+b\hat{p}\right)}$.
We get:

$\displaystyle\exp{(i(X\hat{1}-\mu\hat{q}-\nu\hat{p}))}=\exp{(iX)}\exp{(-i\mu\hat{q})}\exp{(-i\nu\hat{p})}\exp{(i\frac{\mu\nu}{2})}.$
(62)

The matrix element of the latter is

$\displaystyle\Big{[}\exp{(i(X\hat{1}-\mu\hat{q}-\nu\hat{p}))}\Big{]}_{yy^{\prime}}$
$\displaystyle=$
$\displaystyle\exp{(iX)}\exp{(i\frac{\mu\nu}{2})}\Big{[}\exp{(-i\mu\hat{q})}\exp{(-i\nu\hat{p})}\Big{]}_{yy^{\prime}}$
$\displaystyle=$
$\displaystyle\exp{(iX)}\exp{(i\frac{\mu\nu}{2})}\Big{[}\exp{(-i\mu\hat{q})}\Big{]}_{yx}\Big{[}\exp{(-i\nu\hat{p})}\Big{]}_{xy^{\prime}}.$

The matrix element from the latter expression is

$\displaystyle\Big{[}\exp{(-i\mu\hat{q})}\Big{]}_{yx}=\langle
y|\exp{(-i\mu\hat{q})}|x\rangle=\exp{(-i\mu x)}\delta(y-x).$ (64)

The translation operator $\hat{T}(\nu)=\exp{(-i\nu\hat{p}/\hbar)}$ moves
particles and fields by the amount $\nu$, namely

$\displaystyle e^{-i\nu\hat{p}}\Psi(x)=\Psi(x+\nu),\quad
e^{-i\nu\hat{p}}|x\rangle=|x+\nu\rangle.$ (65)

Thus, we have

$\displaystyle\Big{[}\exp{(-i\nu\hat{p})}\Big{]}_{xy^{\prime}}=\langle
x|\exp{(-i\nu\hat{p})}|y^{\prime}\rangle=\langle
x|y^{\prime}+\nu\rangle=\delta(x-y^{\prime}-\nu).$ (66)

Using (5), we can write the matrix element of the density matrix as follows

$\displaystyle\rho(y,y^{\prime})=\frac{1}{2\pi}\int\phi(1;\mu,y-y^{\prime})\exp{\left(-i\frac{\mu(y+y^{\prime})}{2}\right)}d\mu.$
(67)

The conjugate transpose element is

$\displaystyle\rho^{\star}(y^{\prime},y)=\frac{1}{2\pi}\int\phi(-1;\mu,y^{\prime}-y)\exp{\left(-i\frac{\mu(y+y^{\prime})}{2}\right)}d\mu.$
(68)

##  Appendix B Tomogram for the pseudoharmonic oscillator with $a=0$.

Let us first find the tomogram $\mathcal{W}_{n}(X|a,\mu,\nu)$ of (49):

$\displaystyle\mathcal{W}_{n}(X|0,\mu,\nu)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\frac{1}{2\pi|\nu|}\frac{1}{x_{\omega}}\frac{1}{2^{4n+1}}\frac{1}{n!\Gamma(n+\frac{3}{2})}I(X),$
(69) $\displaystyle I(X)$ $\displaystyle\equiv$
$\displaystyle\Bigg{|}\int\limits_{0}^{\infty}\exp{\left[-x^{2}\left(\frac{1}{2x^{2}_{\omega}}-\frac{i\mu}{2\nu}\right)-\frac{iX}{\nu}x\right]}H_{2n+1}\left(\frac{x}{x_{\omega}}\right)dx\Bigg{|}^{2}.$

Using the series expression of the Hermitian polynomial

$\displaystyle
H_{2n+1}\left(\frac{x}{x_{\omega}}\right)=(2n+1)!\sum_{m=0}^{n}\frac{(-1)^{m}}{m!(2n-2m+1)!}\left(\frac{2x}{x_{\omega}}\right)^{2n-2m+1},$
(70)

we can rewrite the integral in the tomogram as

$\displaystyle
I(X)=((2n+1)!)^{2}\Big{|}\sum_{m=0}^{n}\frac{(-1)^{m}}{m!(2n-2m+1)!}\int\limits_{0}^{\infty}e^{\left[-x^{2}\left(\frac{1}{2x^{2}_{\omega}}-\frac{i\mu}{2\nu}\right)-\frac{iX}{\nu}x\right]}\left(\frac{2x}{x_{\omega}}\right)^{2n-2m+1}dx\Big{|}^{2}.$
(71)

We use the integral

$\displaystyle\int\limits_{0}^{\infty}x^{\alpha-1}e^{-px^{2}-qx}dx=\Gamma(\alpha)(2p)^{-\alpha/2}\exp{(q^{2}/8p)}D_{-\alpha}(q/\sqrt{2p}),$
(72)
$\displaystyle[\Re{(\alpha)},\Re{(p)}>0],\quad\text{or}\quad[\Re{(\alpha)},\Re{(q)}>0,\Re{(p)}=0],$
$\displaystyle\text{or}\quad[0<\Re{(\alpha)}<2,\Re{(q)}=\Re{(p)}=0,\Im{(p)}\neq
0],$

where $D_{a}(z)$ is parabolic cylinder functions, to write

$\displaystyle I(X)$ $\displaystyle=$
$\displaystyle((2n+1)!)^{2}e^{-\frac{X^{2}x^{2}_{\omega}}{2\left(\nu^{2}+\mu^{2}x^{4}_{\omega}\right)}}\Bigg{|}\sum_{m=0}^{n}\frac{(-1)^{m}}{m!(2n-2m+1)!}\Gamma(2n-2m+2)$
$\displaystyle\times$ $\displaystyle
D_{-2n+2m-2}\left(\frac{iXx_{\omega}}{\sqrt{\nu}\sqrt{\nu-i\mu
x^{2}_{\omega}}}\right)\left(\frac{2}{x_{\omega}}\right)^{2n-2m+1}\left(\frac{\nu-i\mu
x^{2}_{\omega}}{\nu x^{2}_{\omega}}\right)^{-n+m-1}\Bigg{|}^{2}.$

Then the tomogram is

$\displaystyle\mathcal{W}_{n}(X|0,\mu,\nu)\\!\\!=\\!\\!\frac{1}{\pi}\frac{((2n+1)!)^{2}}{n!\Gamma(n+\frac{3}{2})}\frac{\nu^{2n+1}x_{\omega}}{(\nu^{2}+\mu^{2}x^{4}_{\omega})^{n+1}}e^{-\frac{X^{2}x^{2}_{\omega}}{2\left(\nu^{2}+\mu^{2}x^{4}_{\omega}\right)}}$
$\displaystyle\times$
$\displaystyle\Bigg{|}\sum_{m=0}^{n}\frac{(-1)^{m}}{m!(2n-2m+1)!}\Gamma(2n-2m+2)$
$\displaystyle\times$ $\displaystyle
D_{-(2n-2m+2)}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)\left(\frac{(\nu-i\mu
x^{2}_{\omega})}{4\nu}\right)^{m}\left(\frac{x_{\omega}}{2}\right)^{2m}\Bigg{|}^{2}.$

We use the series representation

$\displaystyle
D_{\nu}(z)=\frac{2^{-\frac{\nu}{2}}e^{-\frac{z^{2}}{4}}}{\Gamma(\frac{1}{2}+\nu)}{}_{1}F_{1}(-\nu,\frac{1}{2},\frac{1}{2}z^{2}),$
(75)

where ${}_{1}F_{1}(a,b,z)$ is the confluent hypergeometric function of the
first kind. It is also known that

$\displaystyle{}_{1}F_{1}(a;b;z)=\sum\limits_{k=0}^{\infty}\frac{(a)_{k}z^{k}}{(b)_{k}k!},$
(76)

where $a_{k}$ is the Pochhammer symbol: $(a)_{k}=a(a+1)\dots(a+k-1)$,
$k=1,2,3\dots$, $(a)_{0}=1$. Combining this knowledge, we can write

$\displaystyle
D_{\nu}(z)=\frac{2^{-\frac{\nu}{2}}e^{-\frac{z^{2}}{4}}}{\Gamma(\frac{1}{2}+\nu)}\sum\limits_{k=0}^{\infty}\frac{(-\nu)_{k}z^{2k}}{{2^{k}}(\frac{1}{2})_{k}k!}=\frac{2^{-\frac{\nu}{2}}e^{-\frac{z^{2}}{4}}}{\Gamma(\frac{1}{2}+\nu)}\sum\limits_{k=0}^{\infty}\frac{\Gamma(-\nu+k)z^{2k}}{\Gamma(-\nu)(2k-1)!!k!},$
(77)

where we used $(x)_{n}=\frac{\Gamma(x+n)}{\Gamma(x)}$ and
$(\frac{1}{2})_{n}=\frac{(2n-1)!!}{2^{n}}$. Then we can write

$\displaystyle D_{-(2n-2m+2)}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)$
$\displaystyle=$
$\displaystyle\frac{2^{n-m+1}e^{\frac{X^{2}x^{2}_{\omega}(\nu+i\mu
x^{2}_{\omega})}{4\nu(\nu^{2}+\mu^{2}x^{4}_{\omega})})}}{\Gamma(\frac{1}{2}-(2n-2m+2))}$
$\displaystyle\times$
$\displaystyle\sum\limits_{k=0}^{\infty}\frac{\Gamma(2n-2m+2+k)\left(\frac{-X^{2}x^{2}_{\omega}(\nu+i\mu
x^{2}_{\omega})}{\nu(\nu^{2}+\mu^{2}x^{4}_{\omega})}\right)^{k}}{\Gamma(2n-2m+2)(2k-1)!!k!},$

holds. It is easy to see that the complex conjugate of the latter function is

$\displaystyle D^{\star}_{-(2n-2m+2)}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)$
$\displaystyle=$ $\displaystyle
D_{-(2n-2m+2)}\left(\frac{iXx_{\omega}\sqrt{\nu-i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right).$ (79)

Then the tomogram is the following

$\displaystyle\mathcal{W}_{n}(X|0,\mu,\nu)\\!\\!=\\!\\!\frac{1}{\pi}\frac{((2n+1)!)^{2}}{n!(n+\frac{1}{2})!}\frac{\nu^{2n+1}x_{\omega}}{(\nu^{2}+\mu^{2}x^{4}_{\omega})^{n+1}}e^{-\frac{X^{2}x^{2}_{\omega}}{2\left(\nu^{2}+\mu^{2}x^{4}_{\omega}\right)}}$
$\displaystyle\times$
$\displaystyle\Bigg{[}\sum_{m_{1},m_{2}=0}^{n}\frac{(-1)^{m_{1}+m_{2}}}{{m_{1}}!{m_{2}}!}\left(\frac{x^{2}_{\omega}}{16\nu}\right)^{{m_{1}}+m_{2}}\left(\nu-i\mu
x^{2}_{\omega}\right)^{m_{1}}\left(\nu+i\mu x^{2}_{\omega}\right)^{m_{2}}$
$\displaystyle\times$ $\displaystyle
D_{-(2n-2{m_{1}}+2)}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)D_{-(2n-2{m_{2}}+2)}\left(\frac{iXx_{\omega}\sqrt{\nu-i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)\Bigg{]}.$

##  Appendix C Tomogram for the pseudoharmonic oscillator for any $a$

Let us find the tomogram for the (5). According to (6) it can be written

$\displaystyle\mathcal{W}(X|\mu,\nu)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi|\nu|}\frac{1}{x_{\omega}}\left[\frac{2n!}{\Gamma(n+\eta+1)}\right]$
$\displaystyle\times$
$\displaystyle\Big{|}\int\limits_{0}^{\infty}\left(\frac{y}{x_{\omega}}\right)^{\eta+\frac{1}{2}}\exp{\left(-\frac{1}{2}\frac{y^{2}}{x^{2}_{\omega}}\right)}L_{n}^{(\eta)}\left(\frac{y^{2}}{x^{2}_{\omega}}\right)\exp{(\frac{i\mu}{2\nu}y^{2}-\frac{iX}{\nu}y)}dy\Big{|}^{2}.$

Changing of variables to $p=\frac{1}{2}-\frac{i\mu}{2\nu}x^{2}_{\omega}$,
$q=\frac{iX}{\nu}x_{\omega}$, we get

$\displaystyle\mathcal{W}_{n}(X|\mu,\nu)=\frac{x_{\omega}}{2\pi|\nu|}\left[\frac{2n!}{\Gamma(n+\eta+1)}\right]\Big{|}\int\limits_{0}^{\infty}t^{\eta+\frac{1}{2}}\exp{\left(-pt^{2}-qt\right)}L_{n}^{(\eta)}(t^{2})dt\Big{|}^{2}.$
(82)

Similarly to the previous subsection, we can use the confluent hypergeometric
function representation:

$\displaystyle
L^{(\eta)}_{n}(t^{2})=\frac{(\eta+1)_{n}}{n}{}_{1}F_{1}(-n;\eta+1;t^{2})=\frac{(\eta+1)_{n}}{n}\sum\limits_{k=0}^{\infty}\frac{(n)_{k}t^{2k}}{(\eta+1)_{k}k!}.$
(83)

Then the tomogram can be written as

$\displaystyle\mathcal{W}_{n}(X|\mu,\nu)=\frac{x_{\omega}}{2\pi|\nu|}\left[\frac{2n!}{\Gamma(n+\eta+1)}\right]\Big{|}\frac{(\eta+1)_{n}}{n}\sum\limits_{k=0}^{\infty}\frac{(n)_{k}}{(\eta+1)_{k}k!}\int\limits_{0}^{\infty}t^{2k+\eta+\frac{1}{2}}\exp{\left(-pt^{2}-qt\right)}dt\Big{|}^{2}.$

We use (72) to write the closed form of the tomogram

$\displaystyle\mathcal{W}_{n}(X|\mu,\nu)=\frac{x_{\omega}}{2\pi|\nu|}\left[\frac{2n!}{\Gamma(n+\eta+1)}\right]\exp{\left(-\frac{X^{2}x^{2}_{\omega}}{2(\nu^{2}+\mu^{2}x^{4}_{\omega})}\right)}\left(\frac{(\eta+1)_{n}}{n}\right)^{2}$
$\displaystyle\times$
$\displaystyle\Bigg{|}\sum\limits_{k=0}^{\infty}\frac{(n)_{k}}{(\eta+1)_{k}k!}\Gamma\left(2k+\eta+\frac{3}{2}\right)\left(\frac{\nu-i\mu
x^{2}_{\omega}}{\nu}\right)^{-(k+\frac{\eta}{2}+\frac{3}{4})}D_{-(2k+\eta+\frac{3}{2})}\left(\frac{iXx_{\omega}\sqrt{\nu+i\mu
x^{2}_{\omega}}}{\sqrt{\nu}\sqrt{\nu^{2}+\mu^{2}x^{4}_{\omega}}}\right)\Bigg{|}^{2}.$


# On wavenumber spectra for sound within subsonic jets

A. Agarwal S. Sinayoko R. D. Sandberg Lecturer, Department of Engineering,
University of Cambridge (aa406@cam.ac.uk)Brunel Research Fellow, ISVR,
University of Southampton (s.sinayoko@soton.ac.uk)Professor, AFM, University
of Southampton<EMAIL_ADDRESS>
###### Abstract

This paper clarifies the nature of sound spectra within subsonic jets. Three
problems, of increasing complexity, are presented. Firstly, a point source is
placed in a two-dimensional plug flow and the sound field is obtained
analytically. Secondly, a point source is embedded in a diverging axisymmetric
jet and the sound field is obtained by solving the linearised Euler equations.
Finally, an analysis of the acoustic waves propagating through a turbulent jet
obtained by direct numerical simulation is presented. In each problem, the
pressure or density field are analysed in the frequency-wavenumber domain. It
is found that acoustic waves can be classified into three main frequency-
dependent groups. A physical justification is provided for this
classification. The main conclusion is that, at low Strouhal numbers, acoustic
waves satisfy the d’Alembertian dispersion relation.

##  1 Introduction

Our initial motivation for understanding the sound spectra in jets came from
the article by Goldstein, (2005) in which he proposed that it may be possible
to identify the “true” sources of noise in jets if the radiating and non-
radiating components could be separated. It is possible to achieve this
separation for the Euler equations linearised about either a steady uniform
base flow (Chu and Kovasznay,, 1958) or a steady parallel flow (Agarwal et
al.,, 2004). Unfortunately, the separation techniques presented in these
papers cannot be applied to full nonlinear Navier-Stokes equations and hence
are not useful for realistic jets.

Sinayoko et al., (2011) showed that filtering in the frequency-wavenumber
domain is an effective technique for separating radiating and non-radiating
components in subsonic jets. Their filtering technique relied on the
dispersion relation $k=|\omega|/c_{\infty}$ (where $k$ denotes the magnitude
of the wavenumber, $\omega$ the angular frequency and $c_{\infty}$ the
farfield speed of sound) satisfied by acoustic waves radiating to a quiescent
farfield. But inside the jet we can have waves that travel supersonically
relative to the ambient medium. In this paper, we define acoustic waves in
jets as those satisfying the dispersion relation
$k_{z}\leq|\omega|/c_{\infty}$, where $k_{z}$ denotes the axial wavenumber. In
other words, in the axial direction, acoustic waves travel either upstream
($k_{z}\leq 0$) or downstream ($0\leq k_{z}\leq|\omega|/c_{\infty}$; in the
downstream case, the axial phase speed is therefore sonic or supersonic. The
characterization of acoustic waves by supersonic axial phase speed was used by
Freund, (2001), Cabana et al., (2008), Tinney and Jordan, (2008) and Obrist,
(2009). The results presented in this paper support this definition of
acoustic waves.

Figure 1: Algorithm for filtering out the radiating field. For numerical
implementation, $W(\bm{k},\omega)$ has a finite width (see Sinayoko et al.,
(2011) for details)

The filtering technique is represented diagrammatically in figure 1. The
radiating part $q^{\prime}(\bm{x},t)$ of a flowfield variable $q(\bm{x},t)$
can be obtained by convolving $q$ with an appropriate filter function
$w(\bm{x},t)$, which is defined in the frequency-wavenumber domain
($W(\bm{k},\omega)$). Sinayoko et al., (2011) considered a model problem in
which the base flow corresponding to the experiment of the Mach 0.9, Re 3600
jet by Stromberg et al., (1980) was excited by two instability waves at
nondimensional frequencies of 2.2 and 3.4. These waves interact nonlinearly to
produce acoustic waves at the difference frequency of 1.2. The density field
at frequency 1.2 is shown in figure 2 (a). In this problem, both acoustic and
hydrodynamic waves are being generated. The Fourier transform of this field,
$P(\bm{k},\omega)$, is shown in figure 2 (b). Multiplying $P(\bm{k},\omega)$
with $W(\bm{k},\omega)$ as defined in figure 1 gives
$P^{\prime}(\bm{k},\omega)$, which is shown in figure 2 (d). The Fourier
transform of the remaining field, $\bar{P}(\bm{k},\omega)=P-P^{\prime}$, is
shown in figure 2 (f). The corresponding density fields in the space-time
domain are obtained by applying the inverse Fourier transforms and are shown
in figures 2 (c) and (e). Radiating components have captured all the acoustic
waves. Clearly the acoustic waves have been separated from the hydrodynamic
waves. However, the efficacy of the filter is puzzling as it is based on the
dispersion relation for sound propagation in a uniform quiescent medium.
Inside the jet we do not have a quiescent medium, so how can this dispersion
relation separate acoustic waves both outside and inside the jet?

|
---|---
|
|

(a) $\rho_{\omega=1.2}(\mathbf{x},t_{0})$ (b) $P(\mathbf{k},1.2)$

(c) $\rho^{\prime}_{\omega=1.2}(\mathbf{x},t_{0})$ (d)
$P^{\prime}(\mathbf{k},1.2)$

(e) $\overline{\rho}_{\omega=1.2}(\mathbf{x},t_{0})$ (f)
$\overline{P}(\mathbf{k},1.2)$

Figure 2: Density field in a turbulent jet of exit radius $a$ plotted in the
physical domain (snapshot in left column) and wavenumber domain (magnitude of
Fourier transform in right column), with colour range $[-5\times
10^{-5},5\times 10^{-5}]$ and $[0,0.5]$ respectively, for a given normalized
frequency $\omega=1.2$ and time $t_{0}$ (c.f. equation (19)). The top row
shows the density field $\rho$, the middle row the radiating field
$\rho^{\prime}$ and the bottom row the non-radiating field $\overline{\rho}$.

In order to answer this question, we have constructed a simple model problem
for sound radiation from a point source in a two-dimensional plug flow (§2).
We show that, for this problem, it is possible to obtain an analytical
expression for the Fourier transform for both the axial and cross-stream
directions. This is a crucial step in obtaining the spectral characteristics
of sound propagation and it enables us to understand and explain the observed
acoustic wavenumber spectra for different frequencies. The solution to this
problem also indicates how to identify acoustic waves for more general
(turbulent) jets. In §3 we consider a more general problem of sound radiation
from a point source in a diverging cylindrical jet and in §4 we identify the
acoustic waves using data obtained from a DNS of a Mach 0.84, Re 7200
turbulent jet.

Even though our motivation for identifying acoustic waves in turbulent jets
stems from a particular application as mentioned above, this work can be used
in other ways. For example, the flow filtering technique defined here could be
used to separate convecting and propagating components in a jet. This can help
define various source models or correlate the nearfield hydrodynamic data to
the acoustic farfield to identify the noise producing regions in the jet. The
technique can also be used to correctly identify the radiating part of
Lighthill’s source term (Freund, (2001), Cabana et al., (2008), Sinayoko and
Agarwal, (2012)).

##  2 Model problem

Figure 3: Schematic sketch of a point source located at the origin inside a
two-dimensional plug flow of width $2a$. The flow Mach number is $M>0$ for
$|y|<a$ (region II) and $M=0$ otherwise (region I).

Consider the problem of a time-harmonic monopole point mass source,
$\rho_{o}\delta(\bm{x})\exp(-i\omega t)$ ($\rho_{o}$ denotes ambient density),
embedded in a plug flow. Several authors (e.g. Morgan, (1975), Mani, (1972))
have considered the problem of farfield sound radiation from a point source in
axisymmetric jets. The main difference between their analysis and ours is that
we seek the spectral content in the frequency-wavenumber domain instead of the
farfield characteristics of sound in the physical domain.

For simplicity, we consider a two-dimensional problem. The problem set up is
described in figure 3. Assuming an $\exp(-i\omega t)$ response
($\phi(\bm{x},t)=\phi(\bm{x};\omega)\exp(-i\omega t)$), the linear velocity
potential $\phi_{I}$ for small disturbances satisfies, in region I (outside
the jet),

$\nabla^{2}\phi_{I}+\kappa^{2}\phi_{I}=0,$ (1)

and in region II, inside the jet,

$\nabla^{2}\phi_{II}-\left(-i\kappa+M\frac{\partial}{\partial
x}\right)^{2}\phi_{II}=\delta(x)\delta(y),$ (2)

where $\nabla^{2}$ denotes the Laplacian operator, $\kappa=\omega/c$ is the
acoustic wavenumber, and $c$ is the speed of sound, which is uniform for the
present problem. Because the velocity potential and pressure are symmetric
about the mid-plane axis of the jet ($y=0$), it is sufficient to solve the
problem for $y\geq 0$. Continuity of pressure at $y=a$ requires that
$p=-\rho_{o}D\phi/Dt$ be continuous (D/Dt denotes material derivative),
therefore

$\left(-i\kappa+M\frac{\partial}{\partial
x}\right)\phi_{II}(x,a)=-i\kappa\phi_{I}(x,a).$ (3)

The kinematic constraint requires that particle displacement $\eta$ at the
interface be continuous. Therefore,

$\frac{\partial}{\partial y}\phi_{I}(x,a)=\frac{\partial}{\partial
t}\eta(x,a),$ (4) $\frac{\partial}{\partial
y}\phi_{II}(x,a)=\frac{\partial}{\partial
t}\eta(x,a)+Mc\frac{\partial}{\partial x}\eta(x,a).$ (5)

Applying the Fourier transform in $x$, defined by

${\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}(k_{x})=\int_{-\infty}^{\infty}\phi(x)e^{-ik_{x}x}dx,$
(6)

to equations (1) and (2), we get,

$\frac{d^{2}{\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{I}}{dy^{2}}+(\kappa^{2}-k_{x}^{2}){\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{I}=0,$ (7) $\frac{d^{2}{\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{II}}{dy^{2}}+\left((\kappa-
k_{x}M)^{2}-k_{x}^{2}\right){\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{II}=\delta(y).$ (8)

Application of the Fourier transform in the axial direction to Eqs. (3) – (5)
gives

$(\kappa-k_{x}M){\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{II}(k_{x},a)=\kappa{\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{I}(k_{x},a),$ (9) $\kappa\frac{d{\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{II}}{dy}(k_{x},a)=(\kappa-
k_{x}M)\frac{d{\hat{\phi\mkern 3.0mu}\mkern-3.0mu}{}_{I}}{dy}(k_{x},a).$ (10)

Let $\beta^{2}=\kappa^{2}-k_{x}^{2}$ and $\gamma^{2}=(\kappa-
k_{x}M)^{2}-k_{x}^{2}$. The locations of the branch cuts for $\beta$ and
$\gamma$ are shown in figure 4. The branch of the square roots are chosen such
that both $\beta$ and $\gamma$ are equal to $\kappa$ for $k_{x}=0$. Acoustic
waves propagate to the farfield only when $\beta$ is real, i.e. when
$|k_{x}|<\kappa$. Therefore, we will focus on this range of wavenumbers. In
region I for outgoing waves to infinity

${\hat{\phi\mkern 3.0mu}\mkern-3.0mu}{}_{I}=Ae^{i\beta y}.$ (11)

Taking into account the symmetry about the mid-plane axis, the solution in
region II is given by

${\hat{\phi\mkern 3.0mu}\mkern-3.0mu}{}_{II}=2B\cos(\gamma
y)-\frac{i}{2}\frac{e^{i\gamma y}}{\gamma}.$ (12)

Application of conditions (9) and (10) yields

$\displaystyle A$ $\displaystyle=$
$\displaystyle-\frac{i\kappa\exp(-ia\beta)(\kappa-
Mk_{x})}{2\Delta(\kappa,k_{x},\gamma)},$ (13) $\displaystyle B$
$\displaystyle=$
$\displaystyle\frac{\exp(ia\gamma)\left(\kappa^{2}\gamma-\beta(\kappa-
Mk_{x})^{2}\right)}{4\gamma i\Delta(\kappa,k_{x},\gamma)},$ (14)

where $\Delta(\kappa,k_{x},k_{y})=\beta\cos(ak_{y})(\kappa-
Mk_{x})^{2}-i\kappa^{2}k_{y}\sin(ak_{y})$. Note that for $M=0$ we recover the
free-field Green’s function of the Helmholtz equation.

Figure 4: Location of the branch cuts for $\beta$ and $\gamma$ (hatched lines)
and the Kelvin-Helmholtz instability pole, $k_{x_{0}}$, in the complex $k_{x}$
domain.

Equation $\Delta(\kappa,k_{x},\gamma)=0$ is the dispersion relation for the
hydrodynamic wave. The roots of this equation represent poles of $\hat{\phi}$
in the complex $k_{x}$ domain. For the present problem there is one root,
$k_{x_{0}}$, which is located in the lower-half $k_{x}$-plane and is
associated with a Kelvin-Helmholtz instability wave; it is purely hydrodynamic
(Agarwal et al.,, 2004) and does not affect our analysis. If our problem had a
pipe (or two splitter plates in 2D) the Kelvin-Helmholtz instability wave
would have an amplitude given by the edge condition, usually an unsteady Kutta
condition (see for example, Crighton, (1985)). Again the Kelvin-Helmholtz wave
would be subsonic and hence would not interfere with our analysis.

Defining the Fourier transform in $y$ by

$\hat{\hat{\phi}}(k_{x},k_{y})=\int_{-\infty}^{\infty}{\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}(k_{x},y)e^{-ik_{y}y}dy,$ (15)

the Fourier transform of the pressure field can be written as

$\displaystyle\hat{\hat{p}}(k_{x},k_{y})=-2i\rho_{o}\int_{0}^{\infty}\left[\kappa{\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{II}(k_{x},y)H(a-y)+\right.$
$\displaystyle\left.(\kappa-k_{x}M){\hat{\phi\mkern
3.0mu}\mkern-3.0mu}{}_{I}(k_{x},y)H(y-a)\right]\cos(k_{y}y)\,dy.$ (16)

This integral can be evaluated analytically:

$\displaystyle\hat{\hat{p}}(k_{x},k_{y})$ $\displaystyle=$
$\displaystyle\frac{\rho_{o}(\kappa-
Mk_{x})}{\Delta(\kappa,k_{x},\gamma)}\left[\frac{i\left(\Delta(\kappa,k_{x},k_{y})-\Delta(\kappa,k_{x},\gamma)\right)}{k_{y}^{2}-\gamma^{2}}+\right.$
(17)
$\displaystyle\left.\frac{\kappa^{2}\left(k_{y}\sin(ak_{y})+i\beta\cos(ak_{y})\right)}{k_{y}^{2}-\beta^{2}}\right].$

###  2.1 Frequency-wavenumber spectra

The acoustic-wave solution in the physical domain can be obtained by applying
the inverse Fourier transforms to Eq. (17) (for details on the geometry of the
Fourier integration contours $F_{x}$ and $F_{y}$ and the implications on
causality, see Agarwal et al., (2004))

$p(x,y;\omega)=\int_{F_{x}}\frac{dk_{x}}{2\pi}e^{ik_{x}x}\int_{F_{y}}\frac{dk_{y}}{2\pi}e^{ik_{y}y}\hat{\hat{p}}(k_{x},k_{y}).$
(18)

If we look at the $k_{y}$ integral, $\hat{\hat{p}}(k_{x},k_{y})$, from Eq.
(17), has two terms. The second term has poles at $k_{y}=\pm\beta$. Using the
method of residues, it can be shown that only these poles contribute to the
integral. Therefore, regardless of the frequency, only the wavenumbers that
satisfy the dispersion relation is $k_{x}^{2}+k_{y}^{2}=\kappa^{2}$ contribute
to the integral. We refer to this as the radiation circle. The other term in
the integrand has zeroes in the denominator at $k_{y}=\pm\gamma$, which
corresponds to an ellipse in the wavenumber domain. Note that these zeroes do
not represent poles as the numerator also goes to zero at $k_{y}=\pm\gamma$.
Therefore, the contribution from the integrand is more complicated for this
term. Further insight can be obtained by plotting the integrand as a function
of frequency. Recall that $\kappa=\omega/c$ and from hereon, for brevity, the
reduced frequency $\kappa a$ is referred to as frequency.

Figure 5 shows the wavenumber spectra $|\hat{\hat{p}}(k_{x},k_{y})|$ for
$M=0.9$ for four different frequencies. At low frequencies ($\kappa a\ll 1$)
most of the energy is concentrated around the radiation circle (figure 5(a)).
For higher frequencies ($\kappa a=O(1)$), the energy is concentrated along the
radiation circle as well, but there is a small amount of energy around the
vertical line $k_{x}=\kappa/(1+M)$ (figures 5(b) and 5(c)). For very high
frequencies ($\kappa a\gg 1$) we see the radiation circle and a part of the
ellipse $k_{y}^{2}=\gamma^{2}$ (figure 5(d)). We observe some ringing around
the ellipse.

For jet noise another useful non-dimensional frequency is the Strouhal number,
$St$, based on the jet diameter and exit velocity. It can be shown that
$St=\kappa a/(\pi M)$. High-speed jet noise peaks at $St\sim 0.2$ ($\kappa
a=0.2M\pi$). This suggests that for filtering out acoustic waves in a jet,
around the peak radiation frequency, one need not worry about the ellipse in
figure 5 (d). For sound radiation near the peak frequency, the dispersion
characteristics are very similar to that of the ordinary wave equation. This
explains why Sinayoko et al., (2011) obtained a good separation of acoustic
and hydrodynamic fields by using a filter based on the dispersion
characteristics of the ordinary wave equation.

#### Low frequencies

A mathematical justification for this low-frequency result can be obtained as
follows. Assume $\kappa a\ll 1$. For acoustic waves, the wavenumbers $k_{x}$
and $k_{y}$ that satisfy the dispersion relation are of the same order of
magnitude as $\kappa$, so that $\gamma a\ll 1$ and $k_{y}a\ll 1$. In Eq. (17)
if we expand the trigonometric functions in a power series up to the second
order, it simplifies to

$\displaystyle\hat{\hat{p}}(k_{x},k_{y})\approx\frac{\rho_{o}\beta(\kappa-
Mk_{x})\left[\kappa^{2}(i-a\beta)+O\left((a\kappa)^{4}\right)\right]}{\Delta(\kappa,k_{x},\gamma)(k_{y}^{2}-\beta^{2})}$

From this equation it is clear that at low frequencies, the dispersion
relation for the problem is $k_{y}^{2}-\beta^{2}=0$, i.e.
$k_{x}^{2}+k_{y}^{2}=(\omega/c)^{2}$, which is the dispersion relation for
sound propagation through a uniform medium at rest. This indicates that mean
flow has a negligible effect on sound propagation at low frequencies. This has
been observed experimentally by Cavalieri et al., (2012).

Figure 5: Wavenumber spectra of the pressure field, $|\hat{\hat{p}}|$, from a
point source in a jet (set-up of figure 3) at frequency, (a) $\kappa
a=0.1\pi$, (b) $\kappa a=\pi$, (c) $\kappa a=2\pi$, (d) $\kappa a=10\pi$. The
density plot uses a linear scale from 0 (white) to $10\pi/\kappa$ (black).

The results can be interpreted by considering the potential field
$\hat{\phi}(k_{x},y)$ given by Eqs. (11) and (12). Figure 6(a) shows the
associated pressure field $\hat{p}(k_{x},y)$ as a function of $y$ for
$k_{x}=\kappa/(1+M)$ for two different values of $\kappa a$, $0.1\pi$ and
$\pi$.

For $k_{x}=\kappa/(1+M)$, $\gamma=0$ and from Eq. (12) it is clear that
$\hat{p}$ does not have a wave-like solution for $y\leq a$. For brevity, we
consider only the real part of $\hat{p}$, which can be shown to be a constant
with respect to $y$ in the limit $\gamma\to 0$. For $y>a$, the response is
given by the harmonic function $\exp(i\beta y)$. At low frequencies (e.g.,
$\kappa a=0.1\pi$), the constant part (when $y<a$) is a very small part of a
wavelength (figure 6(a)) and therefore, when we Fourier transform this field,
most of the energy is concentrated around $k_{y}=\pm\beta$. A similar argument
applies to other values of $k_{x}$. Therefore, the presence of the jet has a
negligible effect on the wavenumber spectrum at low frequencies.

#### Mid-range frequencies

As the frequency increases, the region in the $y$-domain over which the
pressure field is a constant (when $k_{x}=\kappa/(1+M)$) occupies a larger
part of the wavelength (see the dashed line in figure 6 (a)). This has a
significant impact on the Fourier transform (in $y$) of $\hat{p}$. Instead of
being concentrated just around $k_{y}=\pm\beta$, the energy in the $k_{y}$
space gets distributed over a range of $k_{y}$ values. The same reasoning
applies to other values of $k_{x}$ close to $\kappa/(1+M)$. This explains the
vertical patch in the spectra around $k_{x}=\kappa/(1+M)$ in figures 5 (b) and
(c).

Away from these values, in the range $0<k_{x}<\kappa/(1+M)$ , unless the
frequency is very high, the pressure field inside the jet is very similar to
that outside the jet. This is illustrated in figure 6(b), which shows the
pressure field $\hat{p}(k_{x},y)$ for $k_{x}=0.5\kappa/(1+M)$, for $\kappa
a=0.1\pi$ (solid line) and $\kappa a=\pi$ (dashed line). Therefore, for
$0<k_{x}<\kappa/(1+M)$, the Fourier transform of the pressure field is
concentrated around the radiation circle $k_{y}=\pm\beta$.

The energy in the $k_{y}$ direction appears to be contained mainly inside the
radiation circle around $k_{x}=\kappa/(1+M)$. This can be explained as
follows. The difference between the constant field for $y<a$ and the sinusoid
$\exp(i\beta y)$, that would exist in the absence of a jet, is related to the
Heaviside function $H(a-y)$. The Fourier transform of this function is given
by the sinc function, $\mathrm{sinc}(k_{y}a)$. Most of the energy of this
function is contained in the first lobe $k_{y}a<\pi$. This is why there is
little energy outside the radiation circle for $\kappa a=\pi$ (figure 5(b)).
Thus the energy content inside the radiation circle is a consequence of the
pressure field inside the jet. Also the phase speed of the content inside the
radiation circle is supersonic relative to a laboratory reference frame. This
energy content inside the radiation circle represents modes trapped within the
jet. The implication for flow decomposition into radiating and non-radiating
components is that some acoustic components lie within the radiation circle as
well($\kappa^{2}<(\omega/c)^{2}$).

From Eq. (12) and the definition of $\gamma$, it can be seen that for
$k_{x}>\kappa/(1+M)$, we get an exponentially decaying response inside the
jet. Physically this represents a subsonically propagating wave inside the jet
that leads to an evanescent wave. This explains why the spectrum decays for
$k_{x}>\kappa/(1+M)$ for moderate to high frequencies. The decay can also be
explained using ray theory, which predicts a shadow region in a wedge in the
forward direction with half angle $\cos^{-1}[1/(1+M)]$ (Morse and Ingard,
(1968)). A point on the radiation circle determines the direction of sound
radiation to the far field. The angle to the jet axis is given by
$\cos^{-1}(k_{x}/\kappa)$ (Goldstein, (2005)). Thus the shadow region
predicted by ray theory is in agreement with the region of decay on the
radiation circle in Figs. 6(b) – (d).

#### High frequencies

For very high frequencies, we see a combination of the radiation circle and an
ellipse. At high frequencies we would get several wavelengths inside the jet
($y<a$, compare with figure 6(b)). If $a$ is large (say, infinite), then the
spectra would correspond to that of a point source in a uniformly moving
medium, which is an ellipse. This is what we see in figure 5(d). The ringing
is because of the finiteness of $a$, which in the wavenumber domain, results
in the convolution of the sinc function with the radiation ellipse.

Figure 6: Comparison of the pressure field $\hat{p}(k_{x},y)$ at $\kappa
a=0.1\pi$ (solid line) and $\kappa a=\pi$ (dashed line) for (a)
$k_{x}=\kappa/(1+M)$, (b) $k_{x}=0.5\kappa/(1+M)$.

##  3 Diverging jet

In order to consider the effect of three-dimensionality and divergence of a
jet on the results presented in the preceding section, we seek the acoustic
field radiating from a monopole source embedded in diverging axisymmetric mean
flow at $(z/a,r/a)=(3,0)$. The mean flow was obtained from the DNS of a Mach
0.84 and Re 7200 turbulent jet embedded in a co-flow of Mach 0.2. The full
details of the DNS are available in Sandberg et al., (2012), and its sound
field is analysed in the next section. The acoustic field was obtained by
solving the linearised Euler equations.

Figures 7 (a), (c) and (e) show the density field at frequencies of
$ka=0.1\pi$, $\pi$ and $2\pi$, respectively. The notation
$\rho_{ka}(\mathbf{x},t_{0})$ describes a linear combination of the real and
imaginary parts of the temporal Fourier transform of $\rho$ at frequency $ka$,
defined as

$\rho_{ka}(\mathbf{x},t_{0})=\frac{1}{\pi}\left(\Gamma_{r}(\mathbf{x},ka)\cos(\omega
t_{0})-\Gamma_{i}(\mathbf{x},ka)\sin(\omega t_{0})\right),$ (19)

where $\Gamma=\Gamma_{r}+i\Gamma_{i}$ is the temporal Fourier transform of
$\rho(\mathbf{x},t)$. This allows to follow the evolution of the density field
at frequency $ka$ with time. Figures 7 (b), (d) and (f) show the corresponding
wavenumber spectra. Comparing these figures for the case of monopole in the
plug flow at the same frequencies (figures 5(a), (b) and (c)), we can see that
the diverging jet has not changed the nature of the spectra. They look very
similar. We see the radiation circle, the shadow region and the vertical line
around $k_{x}=\kappa/(1+M)$, which can again be identified as trapped waves.
These waves are clearly visible around the centerline in the physical domain
(figure 7 (c) and (e)). The trapped modes propagate to the farfield, as can be
seen by looking at the density along the axis of the jet (figure 8). The decay
rate is a polynomial of order less than 2.5 for all frequencies and not
exponential and thus, these waves are not evanescent.

The presence of trapped modes inside the radiation circle implies that the
filter $|\mathbf{k}|=|\omega|/c_{\infty}$ (figure 1) does not capture all the
acoustic waves within the jet. To do so, one may use the condition
$|\mathbf{k}|\leq|\omega|/c_{\infty}$ instead.

|
---|---
|
|

(a) $\rho_{ka=0.1\pi}(\mathbf{x},t_{0})$ (b) $P(\mathbf{k},0.1\pi)$

(c) $\rho_{ka=\pi}(\mathbf{x},t_{0})$ (d) $P(\mathbf{k},\pi)$

(e) $\rho_{ka=2\pi}(\mathbf{x},t_{0})$ (f) $P(\mathbf{k},2\pi)$

Figure 7: Density field radiating from a monopole in a diverging mean flow of
a turbulent jet plotted in the physical domain (snapshot in left column) and
wavenumber domain (magnitude of Fourier transform in right column), with
colour range $[-5.10^{-5},5.10^{-5}]$ and $[0,0.15]$ respectively, for
normalized frequencies $ka/\pi=0.1,1$ and $2$ and an arbitrary time $t_{0}$
(c.f. equation (19)). In the right column, the arc shows the radiation circle
$|\mathbf{k}|=\omega/c_{\infty}$. The dashed arcs are at
$|\mathbf{k}|=\omega/c_{\infty}\pm\Delta k/2$, where $\Delta k=\sqrt{\Delta
k_{z}^{2}+\Delta k_{r}^{2}}$ is the grid step in the wavenumber domain. Figure
8: Profiles of the density field and the associated decay rate along $r=0$ for
the data of figures 7(a, c, e).

##  4 Turbulent jet

The density field from the DNS Sandberg et al., (2012) of a Mach 0.84 and Re
7200 turbulent jet embedded in a co-flow of Mach 0.2 is shown in figure 9 (a)
at a Strouhal number (based on the jet exit diameter and velocity) of 1.1.
This corresponds to $\kappa a=0.92\pi$.

The usual Fast Fourier Transform algorithm is inconvenient for computing the
temporal Fourier transform, since it requires storing the complete time
history of the three dimensional density field. Here, Goertzel’s algorithm
(1958) was used instead, as it consumes the time history one frame at a time.
No special windowing, i.e. a rectangular window, was used to keep the main
lobe as narrow as possible. However, a small amount of spectral leakage from
low frequencies was observed, resulting in highly supersonic components
($|\mathbf{k}|\leq.6k_{\infty})$ appearing in the wavenumber plots. These
supersonic components are un-physical and have been filtered out by means of a
low pass filter of the form

$W_{L}(\mathbf{k})=\frac{1}{2}\left(1+\tanh((|\mathbf{k}|-k_{0})/(\sigma
k_{\infty})\right),$ (20)

where $k_{\infty}=\omega/c_{\infty}$, $k_{0}/k_{\infty}=0.6$, $\sigma=0.1$.

Figure 9 (b) shows the Fourier transform of the density field. If we filter
only along the radiation circle as for the laminar jet, the density field in
the Fourier domain and the corresponding field in the physical domain are
shown in figures 9 (d) and 9 (c) respectively. A qualitative comparison
between figures 9 (a) and 9 (c) shows that this filter based on the radiation
circle captures most of the acoustic waves. But it can be seen from figure 9
(b) that there is a significant amount of energy inside the radiation circle.
From the above analysis we expect that this region of the spectrum should
contribute to trapped waves along the axis of the jet. Figure 9 (f) shows the
filtered spectra when only the interior of the radiation circle is considered.
The corresponding density field in the physical domain is shown in figure 9
(e). As expected, the waves propagate close to the axis of the jet. These
trapped waves propagate to the far field, as can be seen in the inset of
figure 9 (e), which shows the magnitude of the density field along the
centerline of the jet. As in the preceding section, we get a polynomial decay
rate for these trapped waves.

|
---|---
|
|

(a) $\rho_{St=1.1}(\mathbf{x},t_{0})$ (b) $P(\mathbf{k},1.1)$

(c) $\rho^{\prime}_{St=1.1}(\mathbf{x},t_{0})$ (d)
$P^{\prime}(\mathbf{k},1.1)$

(e) $\overline{\rho}_{St=1.1}(\mathbf{x},t_{0})$

(f) $\overline{P}(\mathbf{k},1.1)$

Figure 9: Density field in a turbulent jet plotted in the physical domain
(snapshot in left column) and wavenumber domain (magnitude of Fourier
transform in right column), with colour range $[-8\times 10^{-6},8\times
10^{-6}]$ and $[0,0.5]$ respectively, for a given normalized frequency
$St=1.1$ and time $t_{0}$ (c.f. equation (19)). The top row shows the density
field $\rho$, the middle row the radiating field $\rho^{\prime}$ and the
bottom row the non-radiating field $\overline{\rho}$. In the right column, the
arc shows the radiation circle $|\mathbf{k}|=\omega/c_{\infty}$. The dashed
arcs are at $|\mathbf{k}|=\omega/c_{\infty}\pm\Delta k/2$, where $\Delta
k=\sqrt{\Delta k_{z}^{2}+\Delta k_{r}^{2}}$ is the grid size in the wavenumber
domain.

Based on the above discussions, we can predict that if we are interested in
the far field acoustic wave radiating at a particular angle $\theta$, then we
can obtain this by filtering on the radiation circle in the frequency-
wavenumber domain for the same polar angle in the $k_{r}-k_{z}$ plane. This is
not surprising and was shown theoretically by Goldstein, (2005). However for
mid to high frequencies, this procedure would miss the effect of trapped waves
that propagate close to the axis of the jet.

It is worth noting that the wavenumber-frequency make-up of the acoustic
spectra is different than that of the turbulence. The turbulent kinetic energy
spectrum is generally obtained from numerical data by the application of the
three-dimensional Fourier transform in space. The resulting spectrum contains
a broad range of wavenumbers and tends to be maximum for low wavenumbers.
Since the radiation circle is also located over low wavenumbers, one may think
that there is no separation between the turbulence and acoustic spectra. This
is not the case because the radiation circle is defined for a particular
frequency: a Fourier transform in time is required in addition to the Fourier
transform in space. In that case, for a given frequency, only a narrow band of
wavenumbers make up the turbulent spectrum. In general, that narrow band would
lie outside the radiation circle (corresponding to subsonic propagation
speeds). For example, to compute the turbulence spectrum from experimental
data, it is common practise to measure at a single point in space and to
Fourier transform the signal in time. The Fourier transform in space is then
obtained by invoking Taylor’s hypothesis: assuming a frozen pattern of
turbulence convected at a local flow speed $U$, the wavenumber and frequency
are related by $\omega/k=U$. Thus, for a given frequency, we are picking out a
single wavenumber of the turbulence spectrum that corresponds to this
dispersion relation. For a subsonic jet this corresponds to a subsonic wave
that lies outside the radiation circle. Thus, for the example problem
considered here the turbulence spectrum corresponds to the energy content to
the right of the radiation circle in 9 (b).

##  5 Conclusions

Most of the acoustic waves radiating from a jet satisfy the d’Alembertian
dispersion relation $k=\omega/c_{\infty}$, i.e. they lie on the radiation
circle in the frequency-wavenumber domain. This validates the radiation
criterion proposed by Goldstein, (2005) and used by Sinayoko et al., (2011).

At low Strouhal numbers (e.g. $St\leq 0.5$), acoustic waves lie mainly on the
radiation circle ($k=\omega/c_{\infty}$). This explains why the dispersion
relation based on an ordinary wave equation was sufficient to filter out the
acoustic waves even inside the jet, shown in figure 2(c).

At mid-Strouhal numbers ($St\sim 1$), some acoustic waves are trapped in the
jet. These trapped waves can be classified as acoustic based on the
observation that:

  1. 1.

they propagate to the far field ;

  2. 2.

they have supersonic phase speed $k<\omega/c_{\infty}$ (in a subsonic jet the
hydrodynamic waves and the energy associated with turbulent structures convect
at subsonic speeds).

These trapped acoustic waves can be identified by using the criterion
$k\leq\omega/c_{\infty}$. Alternatively, one can use the axial wavenumber,
$|k_{x}|\leq\omega/c_{\infty}$ since there are usually no hydrodynamic
components such that $|k_{x}|<\omega/c_{\infty}$ and
$|k_{r}|>\omega/c_{\infty}$. These would represent (unphysical) waves with
axial wavelengths larger than acoustic waves travelling at subsonic speeds at
high angles to the downstream jet axis. An advantage of this approach is that
it does not require computing the radial Fourier transform.

Finally, at high Strouhal numbers ($St\gg 1$), some acoustic waves lie around
the radiation ellipse corresponding to the dispersion relation for waves
propagating through a flow of Mach number equal to the average convection Mach
number. These waves can therefore extend outside the radiation circle.

Although the results presented in this paper are for high Mach number subsonic
jets, the solution of the plug flow problem at low Mach numbers indicates that
the conclusions are valid for low Mach number flows as well.

The above conclusions were shown to hold for sound propagation through a time-
averaged diverging jet and a turbulent jet. Similar results could likely be
obtained for other turbulent flows, such as mixing layers and wakes. If the
flow field in the far field is non-quiescent, which would be the case for a
mixing layer, then the radiation circle turns into a radiation ellipse.

## References

  * Agarwal et al., (2004) Agarwal, A., Morris, P., and Mani, R. (2004).  Calculation of sound propagation in nonuniform flows: suppression of instability waves.  AIAA J., 42(1):80–88.
  * Cabana et al., (2008) Cabana, M., Fortuné, V., and Jordan, P. (2008).  Identifying the radiating core of Lighthill’s source term.  Theoretical and Computational Fluid Dynamics, 22(2):87–106.
  * Cavalieri et al., (2012) Cavalieri, A. V., Jordan, P., Colonius, T., and Gervais, Y. (2012).  Axisymmetric superdirectivity in subsonic jets.  Journal of Fluid Mechanics, 704:388–420.
  * Chu and Kovasznay, (1958) Chu, B. and Kovasznay, L. (1958).  Non-Linear Interactions in a Viscous Heat-Conducting Compressible Gas.  J. Fluid Mech., 3(5):494–514.
  * Crighton, (1985) Crighton, D. G. (1985).  The kutta condition in unsteady flow.  Annual Review of Fluid Mechanics, 17(1):411–445.
  * Freund, (2001) Freund, J. B. (2001).  Noise sources in a low-Reynolds-number turbulent jet at Mach 0.9.  Journal of Fluid Mechanics, 438:277–305.
  * Goertzel, (1958) Goertzel, G. (1958).  An algorithm for the evaluation of finite trigonometric series.  The American Mathematical Monthly, 65(1):34–35.
  * Goldstein, (2005) Goldstein, M. (2005).  On identifying the true sources of aerodynamic sound.  J. Fluid Mech., 526:337–347.
  * Mani, (1972) Mani, R. (1972).  A moving source problem relevant to jet noise.  J. Sound Vib., 25(2):337–347.
  * Morgan, (1975) Morgan, J. (1975).  The interaction of sound with a subsonic cylindrical vortex layer.  Proc. R. Soc. A, 344(1638):341–362.
  * Morse and Ingard, (1968) Morse, P. and Ingard, K. (1968).  Theoretical acoustics.  Princeton University Press.
  * Obrist, (2009) Obrist, D. (2009).  Directivity of acoustic emissions from wave packets to the far field.  Journal of Fluid Mechanics, 640:165–186.
  * Sandberg et al., (2012) Sandberg, R., Suponitsky, V., and Sandham, N. (2012).  DNS of compressible pipe flow exiting into a coflow.  Int. J. Heat Fluid Fl., 35:33–44.
  * Sinayoko and Agarwal, (2012) Sinayoko, S. and Agarwal, A. (2012).  The silent base flow and the sound sources in a laminar jet.  The Journal of the Acoustical Society of America, 131(3):1959.
  * Sinayoko et al., (2011) Sinayoko, S., Agarwal, A., and Hu, Z. (2011).  Flow decomposition and aerodynamic sound generation.  J. Fluid Mech., 668:335–350.
  * Stromberg et al., (1980) Stromberg, J., McLaughlin, D., and Troutt, T. (1980).  Flow field and acoustic properties of a Mach number 0· 9 jet at a low Reynolds number.  Journal of sound and vibration, 72(2):159–176.
  * Tinney and Jordan, (2008) Tinney, C. E. and Jordan, P. (2008).  The near pressure field of co-axial subsonic jets.  Journal of Fluid Mechanics, 611:175–204.