as claimed. ∎

To understand $\epsilon$ in Theorem 13, we must understand the size of
$|\xi|$. To this end, in Section B.3 we will prove that
$|\xi|=\Omega(1.498^{k})$ with constant probability for sufficiently large
$k$. To succeed with probability $1-\delta$ for any $\delta>0$, we note that
the random Gaussian variables appearing in the definition of $|\xi|$ come from
the additional Gaussian variables (the $x^{(\ell)}_{i,j}$ and
$y^{(\ell)}_{i,j}$ terms) employed in the reduction. Therefore, if those
variables are not sufficiently large, we can simply resample them
$O(\log(1/\delta))$-many times until they are.

Let us now turn to our second polynomial interpolation technique, which works
directly with the magnitude of the permanent:

See 14

###### Proof.

We begin by taking the absolute squared value of Equation 47. This leads to an
equation of the form

$\left|\operatorname{Per}(B_{S,T})\right|^{2}=\left|\xi\right|^{2}\left|\operatorname{Per}(A)\right|^{2}+q(z)$
(95)

where $q(z)$ is polynomial in $z$ and $\bar{z}$ with $q(0)=0$. Contrary to
Theorem 13, we cannot compute the left-hand side for $z$ in the unit circle
and leverage the corresponding cancellations, since $q(z)$ has monomials that
depend on $\left|z\right|^{2}$. However, notice that when $z$ is real,
Equation 95 is of the same form as Eq. (8) in [36]. Unsurprisingly, the
remainder of this proof closely follows that in [36], and so, we will just
walk through their proof outline and emphasize the differences that arise in
our case but omit repetition of technical details.

The first step is to choose $x^{(i)}_{i,j}$ and $y^{(i)}_{i,j}$ in $B_{S,T}$
to be i.i.d. Gaussian values from $\mathcal{N}(0,1)_{\mathbb{C}}$ (i.e., the
same distribution as the entries of $A$), and then choose different (real and
positive) values of $z$ in the range $[1-\gamma,1+\gamma]$, for a suitably
small $\gamma$. Note that, if $\gamma=0$, then matrix $B$ has only i.i.d.
elements from $\mathcal{N}(0,1)_{\mathbb{C}}$, as per Theorem 12. When we
choose a different value of $\gamma$ we induce a deformation in the
distribution of $B$, since now some of its elements are drawn from
$\mathcal{N}(0,z^{2})_{\mathbb{C}}$ instead.

If $\gamma$ is too large, then the oracle for
$\left|\operatorname{RGPE}\right|_{\pm}^{2}$ may fail since it only succeeds
with high probability over the Gaussian distribution. To bound how large
$\gamma$ can be before this happen, we compute the total variation distance
between the distributions over $B$ when $\gamma=0$ and when $\gamma>0$ (its
deformed version). If that distance is $O(\delta)$, then the oracle for
$\left|\operatorname{RGPE}\right|_{\pm}^{2}$ will still succeed with
probability $O(\delta)$ and can be invoked to compute
$\left|\operatorname{Per}(B_{S,T})\right|^{2}$ to the desired precision.
Luckily, there is a closed form for the Kullback-Leibler divergence between
two sets of i.i.d. Gaussian variables which, combined with Pinsker’s
inequality, gives a bound on the desired total variation distance (cf. Lemma 1
of [36]). The only thing remaining is to count how many matrix elements of $B$
have actually been deformed via multiplication by $z$. Following the proof of
Theorem 12, we can see that $(c^{2}+k_{S}+k_{T})$ elements of $B$ are standard
complex Gaussians, and the remaining $(c-1)k+k_{S}k_{T}$ have been deformed.
This implies that we must set

$\gamma\leq\frac{\delta}{\sqrt{(c-1)k+k_{S}k_{T}}},$ (96)

for the two distributions to be sufficiently close.

The second step is to apply the polynomial interpolation. This is done by
invoking the oracle to estimate the left-hand side of Equation 95, to error
$\varepsilon(c+k)!\prod s_{i}!t_{j}!$, for $2k+1$ values of $c$ evenly spaced
in the interval $[1-\gamma,1+\gamma]$. This leads to a linear system which can
be solved by a standard least squares method. If we denote by
$\hat{\beta}_{0}$ the resulting estimator for the constant term in the
polynomial of Equation 95, the variance of this estimator is given by

$\operatorname{Var}[\hat{\beta}_{0}]=\frac{\left(\varepsilon(c+k)!\prod
s_{i}!t_{j}!\right)^{2}}{\gamma^{4k}},$ (97)

which is the same as Eq. (20) in [36] when translated to our variables. By
Chebyshev’s inequality, we get

$\Pr\left[\left||\operatorname{Per}(A)|^{2}-\hat{\beta}_{0}|\xi|^{-2}\right|\geq|\xi|^{-2}\sqrt{\operatorname{Var}[\hat{\beta}_{0}]/\delta}\right]\leq\delta,$
(98)

and so we obtain a bound for $\left|\operatorname{Per}(A)\right|^{2}$ that
holds with probability $1-\delta$. We identify this bound with the error in
our estimate for $\left|\operatorname{Per}(A)\right|^{2}$, which needs to be
at most $\epsilon c!$ by the definition of
$\left|\operatorname{GPE}\right|_{\pm}^{2}$. Combining all of the above, we
obtain the following relation

$\epsilon=\varepsilon\frac{(c+k)!\left((c-1)k+k_{S}k_{T}\right)^{k}}{c!X\delta^{2k+1/2}\prod
s_{i}!t_{j}!}\leq\varepsilon\frac{(c+k)^{2k}k^{k}}{X\delta^{2k+1/2}2^{k}},$
(99)

where
$X=\prod_{j=1}^{t_{i}-1}|x^{(i)}_{i,j}|^{2}\prod_{j=1}^{s_{i}-1}|y^{(i)}_{i,j}|^{2}$
is the product of $k$ independent absolute-squared standard i.i.d. complex
Gaussian variables (coming from the definition of $\xi$ in Equation 90). Each
of the $k$ terms, which we will generically call $x$, follows an exponential
distribution with CDF

$\Pr[x\leq g]=1-e^{-g}.$ (100)

If we set $g=\delta/k$, then $\Pr[x\leq g]\leq 1-(1-\delta)^{1/k}$, or in
other words $\Pr[x>g]>(1-\delta)^{1/k}$. By the independence of the $k$ terms
comprising $X$, we get $\Pr[1/X<1/g^{k}]>1-\delta$. Therefore we can assume
that, with probability at least $1-\delta$, $1/X$ is at most
$k^{k}/\delta^{k}$. Plugging this into Equation 99, we get that

$\epsilon=O\biggl{(}\frac{(c+k)^{2k}k^{2k}}{\delta^{3k+1/2}2^{k}}\varepsilon\biggr{)},$
(101)

which is polynomial in $c$, $\delta$ and $\varepsilon$ as long as $k=O(1)$, as
desired. ∎

###  B.3 Lower bound for $|\xi|$

The purpose of this section is to give a lower bound on the $|\xi|$ term
coming from Equation 48 in Theorem 13:

$|\xi|=\prod_{i=1}^{c}s_{i}!t_{i}!\left(\prod_{j=1}^{t_{i}-1}|x^{(i)}_{i,j}|\right)\left(\prod_{j=1}^{s_{i}-1}|y^{(i)}_{i,j}|\right).$
(102)

where $x^{(i)}_{i,j},y^{(i)}_{i,j}$ are i.i.d. standard complex Gaussians. We
show

######  Theorem 26.

$\Pr[|\xi|\geq 1.498^{k}]\geq 1/2-4/\sqrt{k}$.

To prove this theorem, let us first focus on bounding the product of the $k$
Gaussian magnitudes that appear in Eq. (102). To simplify the notation, let us
simply rewrite this product as $X=\prod_{i=1}^{k}X_{i}$ where each $X_{i}$ is
the random variable for the magnitude of a standard complex Gaussian. One can
check that $X_{i}\sim\textrm{Rayleigh}(1/2)$ such that
$\Pr[X_{i}=x]=2xe^{-x^{2}}$ for $x\in[0,\infty)$. To bound $X$, we will need
the following standard central limit theorem:

######  Theorem 27 (Berry-Esseen Theorem).

Let $Z_{1},\ldots,Z_{k}$ be real i.i.d. random variables satisfying

$\displaystyle\mathbb{E}[Z_{i}]$ $\displaystyle=\mu,$
$\displaystyle\mathbb{E}[(Z_{i}-\mu)^{2}]$ $\displaystyle=\sigma^{2}>0,$
$\displaystyle\mathbb{E}[|Z_{i}-\mu|^{3}]$ $\displaystyle=\rho<\infty$

Let $Z=\sum_{i=1}^{k}Z_{i}$ and let $W\sim\mathcal{N}(\mu
k,\sigma^{2}k)_{\mathbb{R}}$ be a real Gaussian. Then, for all
$x\in\mathbb{R}$,

$|\Pr[Z>x]-\Pr[W>x]|\leq\frac{C\rho}{\sigma^{3}\sqrt{k}},$ (103)

where $C$ is some universal constant.

######  Lemma 28.

$\Pr[X\geq.7493^{k}]\geq 1/2-4/\sqrt{k}$.

###### Proof.

Letting $Z_{i}=\log(X_{i})$ and $Z=\sum_{i=1}^{k}Z_{i}$, we begin by writing

$X=e^{\sum_{i=1}^{k}\log(X_{i})}=e^{Z}.$ (104)

One can check that $\Pr[Z_{i}=z]=2e^{2z-e^{2z}}$ leading to the following
moments:

$\displaystyle\mathbb{E}[Z_{i}]$ $\displaystyle=-\gamma/2,$ (105)
$\displaystyle\mathbb{E}[(Z_{i}+\gamma/2)^{2}]$ $\displaystyle=\pi^{2}/24,$
(106) $\displaystyle\mathbb{E}[|Z_{i}+\gamma/2|^{3}]$ $\displaystyle<0.5136$
(107)

where $\gamma\approx 0.5772$ is the Euler-Macheroni constant. Plugging Eqs.
(105), (106), and (107) into the Berry-Esseen Theorem (using $C<.4748$ [59]),
we get that for any $x\in\mathbb{R}$

$|\Pr[Z>x]-\Pr[W>x]|<4/\sqrt{k}$ (108)

where $W$ is real Gaussian with mean $-k\gamma/2$ and variance $k\pi/24$. This
immediately gives the lower bound

$\displaystyle\Pr[Z>-k\gamma/2]$
$\displaystyle=\Pr[Z>-k\gamma/2]-\Pr[W>-k\gamma/2]+\Pr[W>-k\gamma/2]$
$\displaystyle\geq\Pr[W>-k\gamma/2]-|\Pr[Z>-k\gamma/2]-\Pr[W>-k\gamma/2]|$
$\displaystyle\geq 1/2-4/\sqrt{k}$

where the last line uses Eq. (108) and the simple fact that a Gaussian is
greater than its mean with probability $1/2$. Of course, if $Z>-k\gamma/2$,
then $X>(e^{-\gamma/2})^{k}$, which yields the lemma. ∎

Using that $\prod_{i=1}^{c}s_{i}!t_{i}!\geq 2^{k}$ immediately gives Theorem
26.

##  Appendix C Moments of the click distribution

In this appendix we derive the results provided in Section 2.3.3 on click
moments. In Section C.1, we present a brief review of the Gaussian formalism
in terms of covariance matrices an apply it to derive the Husimi covariance
matrix of a Gaussian state in which a transition matrix has been encoded. In
Section C.2, we describe how to reduce the calculation of the mean number of
clicks and its variance for an arbitrary multimode state to the calculation of
one- and two-mode vacuum probabilities. In Section C.3, we show how to fix the
scaling factor $\alpha$ in terms of the photon number density using the
quarter circle law. Finally, in Section C.4, we put together all the results
derived previously and derive Eqs. (24) and (25) giving the first and second
order moments of the click distribution for the two halves of the modes
respectively.

###  C.1 Gaussian states and covariance matrices

In this appendix we provide a self-contained review of the Gaussian formalism
as relevant to the discussion in the main text and the subsequent appendices.
We consider a Gaussian state $\rho$ with $M$ modes and zero mean. The latter
condition means that $\braket{a_{i}}=\text{Tr}(a_{i}\rho)=0$ for all
$i\in\\{1,\ldots,M\\}$. The state $\rho$ is uniquely characterized by its
(complex-)Husimi covariance matrix ${Q}$ with entries

$\displaystyle
Q_{i,j}=\tfrac{1}{2}\braket{z_{i}z_{j}^{\dagger}+z_{j}^{\dagger}z_{i}}+\tfrac{1}{2}\delta_{i,j}=\tfrac{1}{2}\text{Tr}\left(\rho\left[z_{i}z_{j}^{\dagger}+z_{j}^{\dagger}z_{i}\right]\right)+\tfrac{1}{2}\delta_{i,j},$
(109)

where

$\displaystyle{z}=\left(a_{1},\ldots,a_{M},a_{1}^{\dagger},\ldots,a_{M}^{\dagger}\right).$
(110)

For $Q$ to represent a valid quantum mechanical state it needs to satisfy the
uncertainty relation [60]

$\displaystyle Q+\tfrac{1}{2}Z-\tfrac{1}{2}\mathbb{I}_{2M}\succeq 0,$ (111)

where $Z=\left(\begin{smallmatrix}\mathbb{I}_{M}&0\\\
0&-\mathbb{I}_{M}\end{smallmatrix}\right)$, with $\mathbb{I}_{M}$ the $M\times
M$ identity matrix.

The marginal state of a subset of the modes of a Gaussian state is another
Gaussian state. In particular the covariance matrix of a subset of the modes
$i_{1},\ldots,i_{K}$ is obtained by keeping rows and columns
$i_{1},\ldots,i_{K}$ and $i_{1}+M,\ldots,i_{K}+M$ from the original covariance
matrix $Q$. Finally, note that the diagonal elements of the Husimi covariance
matrix are related to the mean photon number in each of the modes as

$\displaystyle Q_{i,i}=Q_{i+M,i+M}=\braket{a^{\dagger}_{i}a_{i}}+1.$ (112)

For a given zero-mean Gaussian state the probability of obtaining vacuum in
all its modes when measuring its photon number statistics is given by

$\displaystyle\Pr(\text{vac})=\frac{1}{\sqrt{\text{det}(Q)}}.$ (113)

A simple expression exists linking a pure Gaussian state parametrized by an
adjacency matrix $A=A^{T}\in\mathbb{C}^{M\times M}$ and its Husimi covariance
matrix [7, 19]

$\displaystyle Q=\left(\mathbb{I}_{2M}-\mathcal{X}(A\oplus
A^{*})\right)^{-1},\text{ where
}\mathcal{X}=\left(\begin{smallmatrix}0&\mathbb{I}_{M}\\\
\mathbb{I}_{M}&0\end{smallmatrix}\right).$ (114)

We now consider the case where $M=2m$ and we use $A$ to encode a transition
matrix $C\in\mathbb{C}^{m\times m}$ as follows

$\displaystyle A=\begin{pmatrix}0&C\\\ C^{T}&0\end{pmatrix}.$ (115)

We write the following decomposition $C=U\sigma V^{T}$ with $U$ and $V$
$m\times m$ unitary matrices and $\sigma$ a diagonal matrix with
$0\leq\sigma_{i,i}<1$. We can then write

$\displaystyle\mathcal{X}(A\oplus A^{*})$
$\displaystyle=\begin{pmatrix}0&0&0&C\\\ 0&0&C^{T}&0\\\ 0&C^{*}&0&0\\\
C^{\dagger}&0&0&0\end{pmatrix}=F\left(\begin{array}[]{cccc}\sigma&0&0&0\\\
0&\sigma&0&0\\\ 0&0&-\sigma&0\\\ 0&0&0&-\sigma\\\
\end{array}\right){F}^{\dagger},$ (120)

where

$\displaystyle{F}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cccc}U&0&-U&0\\\
0&V&0&-V\\\ 0&U^{*}&0&U^{*}\\\ V^{*}&0&V^{*}&0\\\ \end{array}\right),$ (125)

is a unitary matrix of size $4m\times 4m$. With this diagonalization it is
direct to write

$\displaystyle Q=\begin{pmatrix}X&0&0&W\\\ 0&Y&W^{T}&0\\\ 0&W^{*}&X^{*}&0\\\
W^{\dagger}&0&0&Y^{*}\end{pmatrix},$ (126)

with

$\displaystyle
X=U\frac{\mathbb{I}_{m}}{\mathbb{I}_{m}-\sigma^{2}}U^{\dagger}=X^{\dagger},\quad
Y=V\frac{\mathbb{I}_{m}}{\mathbb{I}_{m}-\sigma^{2}}V^{\dagger}=Y^{\dagger},\quad
W=U\frac{\sigma}{\mathbb{I}_{m}-\sigma^{2}}V^{T}.$ (127)

This form of the covariance matrix will be useful below.

###  C.2 Detector click statistics

In this section we derive simple expressions for the mean and variance of the
number of clicks when a quantum state is probed using threshold detectors. We
will obtain expressions for $\braket{c}$, $\braket{c^{2}}$ and $\Delta^{2}c$
that hold for an arbitrary $M$-mode quantum state. In the next subsections we
will specialize these expressions to zero-mean Gaussian states encoding a
random complex-Gaussian transition matrix.

Let $\nu_{c}$ be a subset of $\\{1,2,\dots,M\\}$ where $\nu_{c}$ contains
exactly $c$ elements. Furthermore, let $p(\nu_{c})$ denote the probability of
observing a click in all detectors with labels in $\nu_{c}$. This probability
is given by [61]

$\displaystyle
p(\nu_{c})=\text{Tr}\left[\rho\left(\prod_{j\in\nu_{c}}P_{1}^{(j)}\prod_{j\notin\nu_{c}}P_{0}^{(j)}\right)\right],$
(128)

where $P_{0}^{(j)}=|0_{j}\rangle\langle 0_{j}|$ is a projector onto the vacuum
in mode $j$ and $P_{1}^{(j)}=\operatorname{\mathbbm{1}}^{(j)}-P_{0}^{(j)}$.
The probability of observing $c$ clicks is then given by the sum
$p(c)=\sum_{\nu_{c}}p(\nu_{c})$. Now consider its characteristic function:

$\displaystyle\chi_{c}(x)$
$\displaystyle=\braket{e^{icx}}=\sum_{c=0}^{M}e^{ixc}p(c)$ (129)
$\displaystyle=\sum_{c=0}^{M}e^{ixc}\sum_{\nu_{c}}\text{Tr}\left[\rho\left(\prod_{k\in\nu_{c}}P_{1}^{(k)}\prod_{k\notin\nu_{c}}P_{0}^{(k)}\right)\right]$
(130)
$\displaystyle=\text{Tr}\left[\rho\sum_{c=0}^{M}\sum_{\nu_{c}}\left(\prod_{k\in\nu_{c}}e^{ix}P_{1}^{(k)}\prod_{k\notin\nu_{c}}P_{0}^{(k)}\right)\right]$
(131)
$\displaystyle=\text{Tr}\left[\rho\prod_{\ell=1}^{M}\left(e^{ix}P_{1}^{(\ell)}+P_{0}^{(\ell)}\right)\right]\leavevmode\nobreak\
,$ (132)

where in the last line we used the binomial expansion. The expected number of
detectors that click $\braket{c}$ can then be computed from the characteristic
function [18]

$\displaystyle\braket{c}$
$\displaystyle=-i\frac{d\chi_{c}(x)}{dx}\Bigr{|}_{x=0}$ (133)
$\displaystyle=-i\text{Tr}\left.\left[\rho\frac{d}{dx}\prod_{\ell=1}^{M}\left(e^{ix}P_{1}^{(\ell)}+P_{0}^{(\ell)}\right)\right]\right|_{x=0}$
(134)
$\displaystyle=\text{Tr}\left[\rho\sum_{j=1}^{M}P_{1}^{(j)}\prod_{\ell\neq
j}(P_{1}^{(\ell)}+P_{0}^{(\ell)})\right]$ (135)
$\displaystyle=\sum_{j=1}^{M}\text{Tr}\left[P_{1}^{(j)}\rho\right],$ (136)

where we have used the product rule for derivatives. Using the fact that
$P_{1}^{(j)}+P_{0}^{(j)}=\operatorname{\mathbbm{1}}^{(j)}$, we conclude that

$\braket{c}=M-\sum_{j=1}^{M}\langle P_{0}^{(j)}\rangle,$ (137)

where we have defined $\langle
P_{0}^{(j)}\rangle:=\text{Tr}\left[P_{0}^{(j)}\rho\right]$. Similarly, we can
obtain the second moment of the number of clicks as

$\braket{c^{2}}=\left.(-i)^{2}\frac{d^{2}\chi_{c}(x)}{dx^{2}}\right|_{x=0}.$
(138)

We now write

$\displaystyle\frac{d^{2}\chi_{c}(x)}{dx^{2}}$
$\displaystyle=\sum_{j=1}^{M}\text{Tr}\left(\rho\frac{d}{dx}\prod_{\ell\neq
j}\left[e^{ix}P_{1}^{(\ell)}+P_{0}^{(\ell)}\right]iP_{1}^{(j)}e^{ix}\right)$
(139)
$\displaystyle=\sum_{j=1}^{M}\text{Tr}\left(\rho\Bigr{\\{}i^{2}P_{1}^{(j)}e^{ix}\prod_{\ell\neq
j}\left[e^{ix}P_{1}^{(j)}+P_{0}^{(\ell)}\right]+i^{2}P_{1}^{(j)}e^{ix}\sum_{k\neq
j}^{M}P_{1}^{(k)}e^{ix}\prod_{\ell\neq j,\ell\neq
k}\left[e^{ix}P_{1}^{(\ell)}+P_{0}^{(\ell)}\right]\Bigr{\\}}\right).$

By setting $x=0$ in the last equation we obtain

$\braket{c^{2}}=\sum_{j=1}^{M}\text{Tr}[P_{1}^{(j)}\rho_{j}]+\sum_{j=1}^{M}\sum_{j\neq
k}^{M}\text{Tr}[P_{1}^{(j)}P_{1}^{(k)}\rho].$ (140)

Equivalently, we can use
$P_{1}^{(j)}=\operatorname{\mathbbm{1}}^{(j)}-P_{0}^{(j)}$, which gives

$\braket{c^{2}}=M^{2}-(2M-1)\sum_{j=1}^{M}\braket{P_{0}^{(j)}}+2\sum_{j<k}^{M}\braket{P_{0}^{(j)}P_{0}^{(k)}},$
(141)

where we have defined
$\braket{P_{0}^{(j)}P_{0}^{(k)}}:=\text{Tr}\left[P_{0}^{(j)}P_{0}^{(k)}\rho\right]$.
Putting these results together, the variance of the number of clicks is

$\Delta^{2}c=\sum_{j=1}^{M}\braket{P_{0}^{(j)}}(1-\braket{P_{0}^{(j)}})+2\sum_{j<k}\left(\braket{P_{0}^{(j)}P_{0}^{(k)}}-\braket{P_{0}^{(j)}}\braket{P_{0}^{(k)}}\right)=\sum_{i,j=1}^{M}\left(\braket{P_{0}^{(j)}P_{0}^{(k)}}-\braket{P_{0}^{(j)}}\braket{P_{0}^{(k)}}\right).$
(142)

Since reduced states of Gaussian states can be computed efficiently [62], it
is possible to use these formulas to calculate the average and variance of the
number of clicks for these states. Under the assumption of zero-mean Gaussian
states, the two types of terms in the equations above can be written as

$\displaystyle\braket{P_{0}^{(j)}P_{0}^{(k)}}=\frac{1}{\sqrt{\det(Q^{i,j})}},\quad\braket{P_{0}^{(j)}}=\frac{1}{\sqrt{\det(Q^{j})}},$
(143)

where $Q^{i,j}$ and $Q^{j}$ are the reduced Husimi covariance matrices of
modes $i,j$ and $j$ respectively. We note that these quantities can also be
estimated efficiently for Gaussian states using Monte Carlo methods as
recently shown in Ref. [63].

In the following, we use these results to estimate the average number of
clicks for GBS devices that are configured to encode Gaussian random matrices.

###  C.3 Fixing the scaling factor in terms of the photon number density

In Eq. (14) we used the quarter circle law to find a relation between the
scaling factor $\alpha$, the number of modes $m$ and the total mean photon
number. From this equation we can solve

$\displaystyle\alpha(\mu)=\frac{1+\mu}{\sqrt{\mu}}.$ (144)

For future convenience we can also calculate the following matrix norms of the
blocks of the Husimi covariance matrix (recall Eq. (127)) when the transition
matrix has singular values satisfying the quarter circle law

$\displaystyle\mathbb{E}\left(||X||^{2}\right)$
$\displaystyle=\mathbb{E}\left(||Y||^{2}\right)=\mathbb{E}\left(\sum_{i=1}^{m}\left(\frac{1}{1-\sigma_{i}^{2}}\right)^{2}\right)=m\int_{0}^{2/\alpha(\mu)}p_{\alpha(\mu)}(\sigma)\left(\frac{1}{1-\sigma^{2}}\right)^{2}=m\frac{(1+\mu)}{1-\mu},$
(145) $\displaystyle\mathbb{E}\left(||W||^{2}\right)$
$\displaystyle=\mathbb{E}\left(\sum_{i=1}^{m}\left(\frac{\sigma_{i}}{1-\sigma_{i}^{2}}\right)^{2}\right)=m\int_{0}^{2/\alpha(\mu)}p_{\alpha(\mu)}(\sigma)\left(\frac{\sigma}{1-\sigma^{2}}\right)^{2}=m\frac{\mu(1+\mu)}{1-\mu}.$
(146)

In the equation above we used the fact that the Frobenius norm of a matrix,
$||A||^{2}=\sum_{i,j=1}^{m}|A_{i,j}|^{2}$ can also be expressed as the sum of
the squares of the singular values of the same matrix.

###  C.4 Means and variances of the click distribution

With the preliminary calculations in the last subsections we are ready to
tackle the calculation of the means and covariances of the click distribution
between the two sets of modes when averaged over the set of Gaussian random
matrices.

For this we start with the means, for which we find the covariance of mode $i$
in the first half

$\displaystyle Q^{i}=\begin{pmatrix}X_{i,i}&0\\\
0&X_{i,i}\end{pmatrix}\longrightarrow\braket{P_{0}^{(i)}}=\frac{1}{X_{i,i}}=\frac{1}{1+\braket{a_{i}^{\dagger}a_{i}}},$
(147)

and similarly if $i=j+m$ is in the second half

$\displaystyle Q^{i}=\begin{pmatrix}Y_{j,j}&0\\\
0&Y_{j,j}\end{pmatrix}\longrightarrow\braket{P_{0}^{(i)}}=\frac{1}{Y_{j,j}}=\frac{1}{1+\braket{a_{i}^{\dagger}a_{i}}}.$
(148)

From the last two equations we find
$\braket{d}=m-\sum_{i=1}^{m}\frac{1}{X_{j,j}}$ and similarly
$\braket{e}=m-\sum_{i=1}^{m}\frac{1}{Y_{j,j}}$. Now we can let $X_{j,j}\to
1+\mu$ and $Y_{j,j}\to 1+\mu$ and rearrange to find precisely Eq. (24).

Now we consider the covariances between the clicks in the two halves. For this
we need the two-body covariance matrices of the modes. Consider first the case
where $i$ and $j$ refer to modes in the first half of the set, i.e., $i,j\leq
m$. In that case their reduced covariance matrix is given by

$\displaystyle Q^{i,j}=\begin{pmatrix}X_{i,i}&X_{i,j}&0&0\\\
X_{i,j}^{*}&X_{j,j}&0&0\\\ 0&0&X_{i,i}&X_{i,j}^{*}\\\
0&0&X_{i,j}&X_{j,j}\end{pmatrix}\longrightarrow\braket{P_{0}^{(j)}P_{0}^{(i)}}-\braket{P_{0}^{(j)}}\braket{P_{0}^{(j)}}$
$\displaystyle=\frac{1}{X_{i,i}X_{j,j}-|X_{i,j}|^{2}}-\frac{1}{X_{i,i}X_{j,j}}$
(149)
$\displaystyle=\frac{1}{X_{i,i}X_{j,j}}\sum_{k=1}^{\infty}\left[\frac{|X_{i,j}|^{2}}{X_{i,i}X_{j,j}}\right]^{k}.$
(150)

Note that the Taylor expansion is guaranteed to converge since the uncertainty
relation for the $Q$ covariance matrix (cf. Appendix C.1) guarantees that
$X_{i,i}X_{j,j}>|X_{i,j}|^{2}$. If the two modes $i,j$ are in the second half
one can obtain a similar expression to the one above by letting $X\to Y$

Finally, if $i$ and $j$ are in different halves one has

$\displaystyle Q^{i,j}=\begin{pmatrix}X_{i,i}&0&0&W_{i,j}\\\
0&Y_{j,j}&W_{i,j}&0\\\ 0&W_{i,j}^{*}&X_{i,i}&0\\\
W_{i,j}^{*}&0&0&Y_{j,j}\end{pmatrix}\longrightarrow\braket{P_{0}^{(j)}P_{0}^{(i)}}-\braket{P_{0}^{(j)}}\braket{P_{0}^{(j)}}$
$\displaystyle=\frac{1}{X_{i,i}Y_{j,j}-|W_{i,j}|^{2}}-\frac{1}{X_{i,i}Y_{j,j}}$
(151)
$\displaystyle=\frac{1}{X_{i,i}Y_{j,j}}\sum_{k=1}^{\infty}\left[\frac{|W_{i,j}|^{2}}{X_{i,i}Y_{j,j}}\right]^{k}.$
(152)

As before, the Taylor expansion converges by virtue of the uncertainty
relation.

We can now look at the variance of the total number of clicks (recall Eq.
(142)) over the two halves to write it as

$\displaystyle\Delta^{2}[d+e]=$
$\displaystyle\sum_{i=1}^{m}\left(\frac{1}{X_{i,i}}+\frac{1}{Y_{i,i}}-\frac{1}{X_{i,i}^{2}}-\frac{1}{Y_{i,i}^{2}}\right)$
(153) $\displaystyle+\sum_{i\neq
j}^{m}\left(\frac{1}{X_{i,i}X_{j,j}-|X_{i,j}|^{2}}-\frac{1}{X_{i,i}X_{j,j}}+\frac{1}{Y_{i,i}Y_{j,j}-|Y_{i,j}|^{2}}-\frac{1}{Y_{i,i}Y_{j,j}}\right)$
$\displaystyle+2\sum_{i,j=1}^{m}\left(\frac{1}{X_{i,i}Y_{j,j}-|W_{i,j}|^{2}}-\frac{1}{X_{i,i}Y_{j,j}}\right).$

We will now subtract the quantity
$\sum_{i=1}^{m}\frac{1}{X_{i,i}^{2}}+\frac{1}{Y_{i,i}^{2}}$ in the first term
of the equation above and add it in the second. Then we will also Taylor
expand the terms in the second and third line to write

$\displaystyle\Delta^{2}[d+e]=$
$\displaystyle\sum_{i=1}^{m}\left(\frac{1}{X_{i,i}}+\frac{1}{Y_{i,i}}-\frac{2}{X_{i,i}^{2}}-\frac{2}{Y_{i,i}^{2}}\right)$
(154)
$\displaystyle+\sum_{i,j}^{m}\frac{|X_{i,j}|^{2}}{X_{i,i}^{2}X_{j,j}^{2}}+\frac{|Y_{i,j}|^{2}}{Y_{i,i}^{2}Y_{j,j}^{2}}+\underbrace{\sum_{i\neq
j}^{m}\frac{1}{X_{i,i}X_{j,j}}\sum_{k=2}^{\infty}\left(\frac{|X_{i,j}|^{2}}{X_{i,i}X_{j,j}}\right)^{k}+\sum_{i\neq
j}^{m}\frac{1}{Y_{i,i}Y_{j,j}}\sum_{k=2}^{\infty}\left(\frac{|Y_{i,j}|^{2}}{Y_{i,i}Y_{j,j}}\right)^{k}}_{\equiv\epsilon_{1}}$
$\displaystyle+2\sum_{i,j=1}^{m}\frac{|W_{i,j}|^{2}}{X_{i,i}^{2}Y_{j,j}^{2}}+\underbrace{2\sum_{i,j=1}^{m}\frac{1}{X_{i,i}Y_{j,j}}\sum_{k=2}^{\infty}\left(\frac{|W_{i,j}|^{2}}{X_{i,i}Y_{j,j}}\right)^{k}}_{\equiv\epsilon_{2}}.$

In the last equation we introduced the errors $\epsilon_{1},\epsilon_{2}\geq
0$. We can then write the variance as

$\displaystyle\Delta^{2}[d+e]=$
$\displaystyle\sum_{i=1}^{m}\left(\frac{1}{X_{i,i}}+\frac{1}{Y_{i,i}}-\frac{2}{X_{i,i}^{2}}-\frac{2}{Y_{i,i}^{2}}\right)+\sum_{i,j=1}^{m}\left(\frac{|X_{i,j}|^{2}}{X_{i,i}^{2}X_{j,j}^{2}}+\frac{|Y_{i,j}|^{2}}{Y_{i,i}^{2}Y_{j,j}^{2}}+2\frac{|W_{i,j}|^{2}}{X_{i,i}^{2}Y_{j,j}^{2}}\right)+\epsilon,$
(155)

where $\epsilon=\epsilon_{1}+\epsilon_{2}$. At this point we replace the
diagonal elements $X_{i,i},Y_{j,j}$ by their expectation $1+\mu$ and use the
expressions for the Frobenius norms of $X$, $Y$ and $W$ to find

$\displaystyle\mathbb{E}\left(\Delta^{2}[d+e]\right)\approx
2m\frac{(2-\mu)\mu}{(1-\mu)(1+\mu)^{2}},$ (156)

where we ignore the small error $\mathbb{E}(\epsilon)$.

For the variance in the two halves we can easily write

$\displaystyle\Delta^{2}d=$
$\displaystyle\sum_{i=1}^{m}\left(\frac{1}{X_{i,i}}-\frac{1}{X_{i,i}^{2}}\right)+\sum_{i\neq
j}^{m}\left(\frac{1}{X_{i,i}X_{j,j}-|X_{i,j}|^{2}}-\frac{1}{X_{i,i}X_{j,j}}\right),$
(157) $\displaystyle\Delta^{2}e=$
$\displaystyle\sum_{i=1}^{m}\left(\frac{1}{Y_{i,i}}-\frac{1}{Y_{i,i}^{2}}\right)+\sum_{i\neq
j}^{m}\left(\frac{1}{Y_{i,i}Y_{j,j}-|Y_{i,j}|^{2}}-\frac{1}{Y_{i,i}Y_{j,j}}\right),$
(158)

and use the same approximations used previously to derive the results in Eq.
(25a). Finally, we can obtain the covariance between $d$ and $e$ and the
variance of their difference by recalling that $\Delta^{2}[d\pm
e]=\Delta^{2}d+\Delta^{2}e\pm\text{Cov}(d,e)$.

##  Appendix D Going beyond the dilute limit

The purpose of this section is to give evidence that Theorem 9 can be extended
beyond the dilute limit to show:
$\left|\operatorname{RGPE}\right|_{\pm}^{2}\in\mathsf{FBPP}^{\mathsf{NP}^{\mathcal{O}}}$
for some set of collision patterns that have high probability in the
BipartiteGBS distribution. To do this, we first claim (without
proof)191919This may at first seem counterintuitive since a classical sampler
may choose to favor certain collision patterns over others. However, the
critical point is that the classical sampler does not know which collision
pattern to favor since it will be sampled randomly. that a generalization of
Theorem 7 exists where the subspace $\mathcal{H}_{S}$ can be replaced by any
subspace $\mathcal{H}$ for which $S\in\mathcal{H}$ implies
$S_{\pi,\sigma}\in\mathcal{H}$ for all $\pi,\sigma\in\mathrm{S}_{m}$.
Specifically, let us choose $\mathcal{H}$ to be the space of repetition
patterns with $\mathbb{E}[\braket{n}]$ photon pairs and
$\mathbb{E}[\braket{c}]/2$ clicks in each half of the modes (see Eq. (14) and
Eq. (24), respectively). For this section, we will simply write these
quantities as $n$ and $c$, respectively. Therefore, we have

$|\mathcal{H}|=\binom{m}{c}^{2}\binom{n-1}{n-c}^{2}.$ (159)

Then, the generalization of Corollary 8 allows us to obtain a

$\varepsilon\frac{\mathcal{Z}\alpha^{2n}m^{n}\prod_{i=1}^{m}s_{i}!\prod_{j=1}^{m}t_{j}!}{|{\mathcal{H}}|}$
(160)

approximation to $|\operatorname{Per}(A_{S,T})|^{2}$ where the pattern $(S,T)$
was chosen uniformly from $\mathcal{H}$. On the other hand, to solve
$\left|\operatorname{RGPE}\right|_{\pm}^{2}$, we need an $(\epsilon
n!\prod_{i=1}^{m}s_{i}!\prod_{j=1}^{m}t_{j}!)$-additive approximation.
Following the same logic from Theorem 9, to obtain an
$\mathsf{FBPP}^{\mathsf{NP}^{\mathcal{O}}}$ algorithm, we want

$\mathcal{I}:=\epsilon/\varepsilon=\frac{\mathcal{Z}\alpha^{2n}m^{n}\prod_{i=1}^{m}s_{i}!\prod_{j=1}^{m}t_{j}!}{|{\mathcal{H}}|n!\prod_{i=1}^{m}s_{i}!\prod_{j=1}^{m}t_{j}!}=\frac{\mathcal{Z}\alpha^{2n}m^{n}}{|{\mathcal{H}}|n!}.$
(161)

to be polynomially bounded. Unfortunately, we can no longer use Lemma 5 to
bound this quantity. Because $\alpha=o(m^{1/4})$, the subdominant term of that
expression ($e^{272m/\alpha^{4}}$) would _still_ contribute exponentially to
Eq. (161). Since that expression was obtained through several approximations,
it is no longer tight enough. Nevertheless, we claim that the empirical
scaling of Eq. (161) still appears to be polynomial.

Let us now describe how those numerics were obtained. Given that we do not
have an explicit bound on $\mathcal{Z}$, we must choose to estimate it
somehow. One approach would be to simply draw a random Gaussian matrix and
compute $\mathcal{Z}$ from the eigenvalues of that matrix. Unfortunately, this
approach is quite slow in practice and we would only be able to calculate
$\mathcal{Z}$ for relatively modest numbers of modes. Instead, we will choose
to estimate $\mathcal{Z}$ by an explicit formula. First, consider the quantity

$\mathcal{Z}^{-1}=\det(\mathbb{I}_{m}-A)=\prod_{i=1}^{m}(1-\lambda_{i})$ (162)

for complex Wishart matrix $A\sim\mathcal{W}$ (see Appendix A for explanation
of this ensemble). In other words, the _inverse_ of $\mathcal{Z}$ is the
characteristic polynomial of $A$. The expectation of the characteristic
polynomial for a Wishart matrix is well-understood:

######  Lemma 29 (Edelman [64]).

$\mathbb{E}[\mathcal{Z}^{-1}]=m!\sum_{i=0}^{m}\binom{m}{i}\frac{(-1)^{i}}{(m-i)!(\alpha^{2}m)^{i}}$.

In Figure 4, we show that the quantity $1/\mathbb{E}[\mathcal{Z}^{-1}]$ serves
as a good estimate for $\mathcal{Z}$ by comparing it against its true value
(computed directly with the eigenvalues). In Figure 5, we compute
$\mathcal{I}$ from Eq. (161) using the estimate
$1/\mathbb{E}[\mathcal{Z}^{-1}]$ for $\mathcal{Z}$ and find a relatively
benign scaling. We conclude by reiterating that if we could _prove_ the
accuracy of our estimator and that this scaling is only polynomial, then we
arrive at our desired conclusion:
$\left|\operatorname{RGPE}\right|_{\pm}^{2}\in\mathsf{FBPP}^{\mathsf{NP}^{\mathcal{O}}}$.
Therefore, if one were to also assume that
$\left|\operatorname{RGPE}\right|_{\pm}^{2}$ is $\\#\mathsf{P}$-hard, we
obtain approximate average-case hardness for GBS in the high-collision
$m=\Theta(\mathbb{E}[\braket{n}])$ regime.

$0$$10$$20$$30$$40$$50$$60$$70$$0$$50$$100$$150$$200$$250$$300$$350$$400$$450$$500$

$\log(\mathcal{Z})$

$m$$\alpha=3$$\alpha=2m^{1/8}$$\alpha=2m^{1/4}$ Figure 4: Comparison of
$1/\mathbb{E}[\mathcal{Z}^{-1}]$ and $\mathcal{Z}$ as a function of $m$ and
various $\alpha$ (shown on a log scale). Each point was calculated from Lemma
29, and the error bars show the sample standard deviation of
$\log(\mathcal{Z})$ calculated from 50 complex Wishart matrices drawn from
$\mathcal{W}(0,\frac{1}{\alpha^{2}m}\mathbb{I}_{m})$.

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

$\mathcal{I}$(millions)

$m$(thousands)$\alpha=3$$\alpha=4$

$9$$10$$11$$12$$13$$14$$15$$16$$17$$9.5$$10$$10.5$$11$$11.5$

$\log(\mathcal{I})$

$\log(m)$$\alpha=3$(slope=1.500017)$\alpha=4$(slope=1.500036)

Figure 5: Left:
$\mathcal{I}:=\frac{\mathcal{Z}\alpha^{2n}m^{n}}{|{\mathcal{H}}|n!}$ as a
function of $m$ and various $\alpha$ (domain and range have been scaled by
$10^{3}$ and $10^{6}$, respectively). Large values of $\alpha$, such as those
corresponding to the dilute limit are dwarfed by the scaling of $\mathcal{I}$
at these small values of $\alpha$, and so are not shown. Right: The same data
as the left graph, except shown on a log-log plot and only values of $m$
greater than $10^{4}$. Linear regression lines suggest that $\mathcal{I}$
scales like $m^{3/2}$ for small constant $\alpha$.

###  D.1 Justification for collision subspace

The purpose of this section is to justify that the subspace $\mathcal{H}$ only
contains collision patterns that suffice for the reduction given by Lemma 24.
While it is known that Lemma 24 is insufficient to prove an efficient
reduction between $\operatorname{GPE}_{\pm}$ and $\operatorname{RGPE}_{\pm}$,
it is still weak evidence that the permanents that appear in the probabilities
are hard to estimate. To recap that procedure, we embed the permanent of
matrix $A$ into the permanent of a matrix $B_{S,T}$ with repeated rows/columns
provided the repetition patterns $S$ and $T$ have a least $k_{T}$ and $k_{S}$
modes (respectively) with a single photon.202020Recall that while $k_{T}$ and
$k_{S}$ are the number of collisions, they are also the number of modes that
we need to have without collisions. This is due to the fact that every
repeated row/column requires an _unrepeated_ column/row in the reduction.
Notice that $\mathcal{H}$ contains repetition patterns with $k_{T}=k_{S}=n-c$.
So, to match the reduction, $n-c$ of the $c$ modes that click must have
exactly one photon. We will show that under a relatively modest condition on
$\alpha$, all repetition patterns in $\mathcal{H}$ will satisfy this property.

While this may seem like a rather strong condition, the fact that the number
of modes that click is always a significant fraction of the total number of
photons makes it possible. First notice that there are at most $n-c$ modes (of
our total budget of $c$ clicks) that can have more than two photons. For the
reduction, we then need a different set of $n-c$ modes to click. That is, we
want

$2(n-c)\leq c.$ (163)

Using Eq. (14) and Eq. (24), this will hold whenever $\alpha\geq 3/\sqrt{2}$.

Finally, we claim that $A$ will be a $(2c-n)\times(2c-n)$ matrix because the
total number of clicks must be the dimension of $A$ plus the $n-c$ extra
clicks required for the reduction (i.e., $c=(2c-n)+(n-c)$). In order to obtain
a hardness proof, we would like $A$ to be relatively large so that computing
its permanent is $\\#\mathsf{P}$-hard. We have that the size of $A$ is at
least half the total number of clicks in the distribution whenever $2c-n\geq
c/2$. Notice that this is the same condition as that in Eq. (163), and so we
only require $\alpha\geq 3/\sqrt{2}$. This is more than sufficient to
guarantee hardness for all regimes from the dilute limit all the way to when
the number of photons is linear in the number of modes. We stress again,
however, that our reduction is only efficient given an oracle for
$\operatorname{Per}(B_{S,T})$, and the numerics of the previous section only
suggest an oracle for $|\operatorname{Per}(B_{S,T})|^{2}$.


# Caustics type of gravitational singularities

G. A. Sardanashvily  E. G. Timoshenko Translated on 07.07.2024 by
<EMAIL_ADDRESS>from Russian: The Bulletin of Moscow State
University, ISSN 2074-6636, UDC 530.12, Series 3, Physics and astronomy, Vol.
30, No 3, pp. 75-77 (1989).

( Division of Theoretical Physics, Department of Physics, Lomonosov Moscow
State University, Russia, USSR )

###### Abstract

We describe a new type of gravitational singularities which are caustics of
spatial–temporal foliations.

##  1 Introduction

The description of gravitational singularities as singularities of space–time
foliations [1-3] allows one, unlike other methods, to establish the structure
of these singularities. In this paper, we consider the caustics of
spatial–temporal foliations.

There is a diagram

$\begin{array}[]{ccc}GL(4,\mathbb{R})&\longrightarrow&SO(4)\\\
\downarrow&&\downarrow\\\ SO(3,1)&\longrightarrow&SO(3)\end{array}$

of the reduction of the structural group of the tangent bundle $T(X^{4})$ over
a smooth orientable paracompact manifold $X^{4}$ admitting the gravitational
field $g$ as a global section of the bundle of pseudo–Euclidean bilinear forms
in tangent spaces to $X^{4}$. By virtue of this diagram, any gravitational
field on the manifold $X^{4}$ defines some smooth orientable space–time
distribution $F$ of 3-dimensional space–like (with respect to $g$) hyperplanes
of tangent spaces to $X^{4}$, given by the equation $h^{0}(F)=0$, where
$h^{0}=h^{0}_{\mu}dx^{\mu}$ is the tetrad form of the field $g$. Inversely,
any 3-dimensional smooth orientable distribution on $X^{4}$ is a space–time
with respect to some gravitational field [2,3].

A distribution $F$ is called integrable (foliation) if its generating form
$h^{0}$ satisfies the equation $h^{0}\wedge dh^{0}=0$. In this case, the
manifold $X^{4}$ splits into 3-dimensional space–like hypersurfaces such that
the layers of the distribution are tangent spaces to these hypersurfaces. A
space–time foliation is said to be causal if $h^{0}=Ndf$, where $N$ and $f$
are real functions on $X^{4}$ such that $N\not=0$ and $df\not=0$ are
everywhere on $X^{4}$.

The correspondence between gravitational fields and spacetime distributions on
the manifold $X^{4}$ allows us to describe gravitational singularities as
singularities of such distributions. Since any regular gravitational field
locally admits a causal foliation, singularities of distributions can be
described as singularities of causal foliations in terms of generating
functions $f$. Two types of such singularities are distinguished.

1\. The generating function $f$ is single–valued but has a critical point
where $df=0$. It gives rise to a Hefliger structure (singular foliation of $f$
level surfaces on $X^{4}$, which change their topology when passing through
the critical point of the function $f$.

2\. The function $f$ is multivalued. In the region where it is single-valued,
the function $f$ generates a foliation, which, however, collapses at branching
points 1, where the layers begin to intersect each other. To describe this, we
can lift the foliation $F$ into the total space of the cotangent bundle
$\mbox{tl}\ T^{*}(X^{4})$, continue it there and then project it back to the
base $X^{4}$. The singularities of $F$ will correspond to the singularities of
this projection, which are caustic points.

In gravitation theory (by analogy with geometrical optics), a caustic is
usually a geometric place of focal or conjugate points [4,5]. We follow a more
general definition of caustics as singularities of Lagrangian mappings [6].
Any such caustics can be locally (in terms of germs) reduced to the following
canonical construction.

Consider space $R^{2n}$ with coordinates $(x^{\mu},P_{\mu})$ 1- form
$a=P_{\mu}dx^{\mu}$ on $R^{2n}$ and a $n$-dimensional submanifold $N\subset
R^{2n}$ such that $da|_{N}=0$. This is called Lagrangian submanifold and it
can be defined by the generating function $S(x^{i},P^{j})$ of $n$ variables
$(x^{i},P_{j},\ i\in I,j\in J)$, where $(I,J)$ is some partition of the set
($1$,…$n$), and is defined by the relations

$x^{j}=-\frac{\partial S}{\partial P_{j}},\quad P_{i}=\frac{\partial
S}{\partial x^{i}}.$

Let $\pi:\ \\{x^{\mu},P_{\mu}\\}\to\\{x^{\mu}\\}$ be the projection of
$R^{2n}$ onto $R^{n}$. Being restricted on the Lagrangian manifold:

$\pi_{N}:(x^{i},P_{j})\rightarrow(x^{i},x^{j}=-\partial S/\partial P_{j}),$

it is called a Lagrangian mapping, and the set of its special points (where
the matrix $\partial^{2}S/\partial P_{i}\partial P_{j}$ is degenerate) is
called a caustic.

Let $f$ be the generating function of the foliation $F$ on $X$. Let us define
the embedding

$\gamma:\\{x^{\mu}\\}\to\\{x^{\mu},P_{\mu}=\partial f/\partial x^{\mu}\\}$

of manifold $X$ into the total space tl$\,T^{*}(X)$. Its image is a Lagrangian
submanifold of tl$\,T^{*}(X)$ endowed with an induced foliation
$F^{\prime}=\pi^{*}_{\gamma(X)}F$, where $\pi_{\gamma(X)}:\ \gamma(X)\to X$ is
a Lagrangian mapping. Inversely, let $N\subset\mbox{tl}\,T^{*}(X)$ be a
Lagrangian submanifold with a generating function $S(x^{i},P_{j})$. Let
$F^{\prime}$ denote the foliation of level surfaces of the function

$f^{\prime}(x^{i},P_{j})=S-P_{j}\frac{\partial S}{\partial P_{j}}$

on $N$. The image $\pi_{N}(F^{\prime})$ of the foliation $F$ under the
Lagrangian mapping $\pi_{N}$ forms some structure on $X$, which is a foliation
with the generating function $f(x^{\mu})=f^{\prime}(x^{i},P_{j}(x^{i},x^{j}))$
on the region where $\pi_{N}$ has no singularities. This foliation breaks down
at the caustic points of the mapping $\pi_{N}$, where the functions
$P_{j}(x^{i},x^{j})$ and $f(x^{\mu})$ become multivalued. The classification
of stable caustics $A_{2-5},\ D_{4,5}$ is known on a 4-dimensional manifold.

Gravitational singularities of the caustics type have the following
distinctive feature. There exist regions where not the nearest, but, on the
contrary, the layers of foliation far separated by the values of $f$
intersect. Therefore, the spatial–temporal foliation and, consequently, the
gravitational field in such a region can be locally extended beyond the
caustic points. Let us give the following example.

We define $f(u,v)$ as a function on $R^{2}$ satisfying the equation
$f^{3}-3uf-2v=0$. This function is single-valued:

$f_{+}=[v+(v^{2}-u^{3})^{1/2}]^{1/3}+[v-(v^{2}-u^{3})^{1/2}]^{1/3}$

in the region $U=(u,v:v^{2}>v^{3})$ and is three–valued:

$f_{0,1,2}=2u^{1/2}\cos\left(\frac{1}{3}(\varphi+2\pi
n)\right),\quad\varphi=\arccos(vu^{-3/2}),\quad n=0,1,2$

at points $v^{2}<u^{3}$. Let $F$ be a foliation $f_{+}(u,v)=c=const$ on
$U\subset R^{2}$. Its layers $F_{c}$ are straight lines $2v=c^{3}-3cu$. This
foliation forms a caustic at the branching points $v^{2}=u^{3}$, so that
$u=v=0$ is a caustic of type $A_{3}$ whereas $v^{2}=u^{3}\not=0$ $A_{2}$-type
caustic points. The layers $F_{c}$, $c>\varepsilon>0$ can be extended beyond
the caustic line $v=u^{3/2}$ as $f_{0}(u,v)=const$ layers in the region
$0<v<u^{3/2}$. But they begin to intersect each other when $v<0$, although the
layers closest to each other intersect only near the caustic line
$v=-u^{3/2}$. At the same time, the layers $F_{c>0}$ begin to intersect with
the layers $F_{c<0}$ already at the caustic line $v=u^{3/2}$.

An example of singularity with caustics of the above type gives the
gravitational wave $(h^{0}=Nd(f_{+}+x^{2}))$

$ds^{2}=2dx^{2}du+\left[-f_{+,v}^{2}-(f_{+,u}-1)^{2}\right]du^{2}-dv^{2}-(dx^{3})^{2}.$

## References

  * [1] D. Ivanenko, G. Sardanashvily // Phys. Lett. 1982. 91A. P. 341.
  * [2] G. Sardanashvily // Acta Phys. Hung. 1985. 57. P. 31.
  * [3] G. Sardanashvily, V. Yanchevsky // Acta Phys. Pol. 1986. B17. P. 1017.
  * [4] F. Warner // Am. J. Math. 1965. 87. P. 575.
  * [5] K. Rosquist // Comm. Math. Phys. 1983. 88. P. 339.
  * [6] V. I. Arnold, S. M. Gusein-Zade, A. N. Varchenko // Singularities of Differentiable Maps. 1988. Birkhäuser, Boston. ISBN 978-1-4612-8408-6.

Manuscript (Brief Communications, Theoretical and Mathematical Physics)
received in the editorial office on: 23.06.88.