# Quantum dynamics of gravitational massive shell

Andrzej Góźdź<EMAIL_ADDRESS>Institute of Physics, Maria Curie-
Skłodowska University, pl. Marii Curie-Skłodowskiej 1, 20-031 Lublin, Poland
Marcin Kisielowski<EMAIL_ADDRESS>Department of Fundamental
Research, National Centre for Nuclear Research, Pasteura 7, 02-093 Warszawa,
Poland  Włodzimierz Piechocki<EMAIL_ADDRESS>Department of
Fundamental Research, National Centre for Nuclear Research, Pasteura 7, 02-093
Warszawa, Poland

(August 28, 2024)

###### Abstract

The quantum dynamics of a self-gravitating thin matter shell in vacuum has
been considered. Quantum Hamiltonian of the system is positive definite.
Within chosen set of parameters, the quantum shell bounces above the horizon.
Considered quantum system does not collapse to the gravitational singularity
of the corresponding classical system.

###### Contents

  1. I Introduction
  2. II Classical dynamics
  3. III Quantum level
    1. III.1 Affine coherent states
    2. III.2 Quantum observables
    3. III.3 Quantum dynamics
  4. IV Matrix elements of quantum observables
    1. IV.1 New expression for the basis elements
    2. IV.2 The matrix elements of the Hamiltonian operator
    3. IV.3 The matrix elements of the momentum operator
      1. IV.3.1 Entries below and on the diagonal
      2. IV.3.2 The entries above the diagonal
      3. IV.3.3 Summary
    4. IV.4 The matrix elements of the coordinate operator
      1. IV.4.1 Diagonal terms
      2. IV.4.2 The terms above the diagonal
      3. IV.4.3 Summary
  5. V Evolution of quantum observables
    1. V.1 Gaussian state
    2. V.2 Eigenvalues convergence
    3. V.3 Evolution of the observables
    4. V.4 Oscillatory character of the evolution
  6. VI Conclusions
  7. A Orthonormal basis of the carrier space

##  I Introduction

The dynamics of a self-gravitating thin matter shell is one of the simplest
models describing gravitational collapse of an isolated gravitational system.
In the case of a spherically symmetric shell in vacuum, a satisfactory
Hamiltonian description of that dynamics has been found; see the paper [1] and
references therein. By shell in vacuum one means a thin matter shell with a
region of flat Minkowski space in the interior and the Schwarzschild geometry
in the exterior of the shell. The global Hamiltonian of that system is
explicitly time independent and is a function of two canonically conjugated
phase space variables. That Hamiltonian is equal to the Arnowit-Deser-Misner
(ADM) mass at spacial infinity. Having well defined Hamiltonian description of
matter shell, we have decided to quantize that system to get insight into
corresponding quantum dynamics. Present paper is devoted to the examination of
such an issue.

The shell system is simple enough to be treated satisfactory at classical
level, and rich enough for the examination of various aspects of corresponding
quantum system. Recently, it was used for addressing the issue of the
importance of the choice of time parameter at quantum level [2]. It was shown
that quantum theories of the shell for different choices of time are not
unitarily equivalent.

We have found that our quantum Hamiltonian is positive definite which supports
its classical property. Within chosen set of parameters describing our system,
the quantum shell bounces above the horizon that is in contrast to the
classical case. It means that due to quantum effects our quantum system does
not collapse to gravitational singularity.

The paper is organized as follows: In Sec. II we present the solution to
Hamilton’s dynamics restricting considerations to the subspace of phase space
for which the Hamiltonian is positive definite. In Sec. III we recall the
coherent states quantization method applied in this paper. Sec. IV concerns
the calculations of the matrix elements of quantum observables in specific
basis of considered Hilbert space. That includes the operator of Hamiltonian
and operators of canonical variables. The quantum evolution of the system is
presented in Sec. V. We conclude in Sec. VI.

In the following we choose $\;G=c=1=\hbar\;$ except where otherwise noted.

##  II Classical dynamics

For self-consistency of the present paper, we recall the main results of Ref.
[1]. Next, we present the solution to the classical dynamics.

The canonical structure of the phase space of the system “shell+gravity” is
given by

$\omega=dp\wedge dq\,,$ (1)

where $q\in\mathbb{R}_{+}:=\\{x\in\mathbb{R}\;|\;x>0\\}$ is the configuration
variable representing the proper volume of the shell, and $p\in\mathbb{R}$ is
the momentum representing the hyperbolic angle between the surfaces of
constant time on both sides of the shell.

Hamilton’s dynamics reads

$\displaystyle\dot{q}=\frac{\partial H}{\partial p}\,,$ (2)
$\displaystyle\dot{p}=-\frac{\partial H}{\partial q}\,,$ (3)

where the Hamiltonian is defined to be

$H(p,q):=\sqrt{\frac{q}{2}}\;\Big{[}1-\Big{(}\cosh(p)-\sqrt{\frac{m^{2}(q)}{2q}+\sinh^{2}(p)}\Big{)}^{2}\Big{]}\,,$
(4)

and where $m(q)$ represents the total rest mass of the matter (energy) of the
shell and plays the role of the constitutive equation for the matter field of
the shell. The dots over $q$ and $p$ in (2)–(3) denote time derivatives, where
the time variable is the Schwarzschild time $t$ measured at spatial infinity.
Since $H$ is time independent, the total energy of the entire system is
conserved.

Making use of (2)–(4) we can determine if the shell’s size increases or
decreases with time. It follows that

$\dot{q}=\frac{\partial H}{\partial
p}=\sinh(p)\sqrt{2q}\frac{(\cosh(p)-\sqrt{\frac{m^{2}}{2q}+\sinh^{2}(p)})^{2}}{\sqrt{\frac{m^{2}}{2q}+\sinh^{2}(p)}}.$
(5)

Therefore the sign of $\dot{q}$ is dictated (up to singularities where the
right-hand-side vanishes) by $\sinh(p)$. This means that for positive $p$ the
shell grows $\dot{q}\geq 0$ and for negative $p$ the shell shrinks
$\dot{q}\leq 0$.

Any canonical transformation of the system (1)–(4) leads to the physically
equivalent system. The advantage of the present choice of the phase space
variables is their clear physical interpretation.

The Hamiltonian function (4) is equal to the total energy of considered
isolated gravitational system at spatial infinity so that it is the ADM mass.
It means, roughly speaking, that if the density of considered matter field is
positive, the ADM mass must be positive (see, [3] and references therein).
This may impose the restriction on the phase space of considered gravitational
system. The specific form of matter field may lead to the specific subspace
$\Lambda$ of the phase space $\Pi=\\{(p,q)\leavevmode\nobreak\
|\leavevmode\nobreak\ p\in\mathbb{R},q\in\mathbb{R}_{+}\\}$, such that
$H(p,q)>0$ for $(p,q)\in\Lambda$. We call $\Lambda$ the physical phase space.

Let us consider the case $m(q):=m=const$ which corresponds to the dust matter.
Since the Hamiltonian is time independent, the energy is conserved. For each
value of the Schwarzschild mass $M>0$, the curve

$H(p,q)=M$ (6)

is the shell’s trajectory in the phase-space. Let us notice, that the
Hamiltonian is an even function of $p$ and therefore it is enough to look for
solutions with positive momentum $p>0$. Rewriting the equation (6) (see also
(4)) in the form

$1-\sqrt{\frac{2}{q}}M=\;\Big{(}\cosh(p)-\sqrt{\frac{m^{2}(q)}{2q}+\sinh^{2}(p)}\Big{)}^{2},$
(7)

we immediately notice that $q>2M^{2}$. This condition says that the proper
volume $q$ is bounded from below so that reflects the fact that the shell is
outside the event horizon of the exterior Schwarzschild solution.

The equation (6) can be solved for $p$ as a function of $q$ and $M$. It can be
checked by substitution that

  * •

for $M\sqrt{\frac{2}{q}}\geq\frac{m^{2}}{2q}$, i.e.,
$\sqrt{q}\geq\frac{m^{2}}{2\sqrt{2}M}$, the solution satisfies:

$\cosh(p)=\frac{2-\frac{m^{2}}{2q}-M\sqrt{\frac{2}{q}}}{2\sqrt{1-M\sqrt{\frac{2}{q}}}}=:f(q),$
(8)

  * •

while for $M\sqrt{\frac{2}{q}}<\frac{m^{2}}{2q}$,
i.e.,$\sqrt{q}<\frac{m^{2}}{2\sqrt{2}M}$, it satisfies:

$\cosh(p)=-\frac{2-\frac{m^{2}}{2q}-M\sqrt{\frac{2}{q}}}{2\sqrt{1-M\sqrt{\frac{2}{q}}}}=-f(q).$
(9)

At the first glance, it may seem that for $0\leq M<\frac{m}{2}$ both relations
(8) and (9) are relevant and that there is a discontinuity at
$q=\frac{m^{2}}{2\sqrt{2}M}$. However, it turns out that in this case for
$\sqrt{q}\geq\frac{m^{2}}{2\sqrt{2}M}$ the function $f(q)$ takes only values
smaller than $1$ and therefore the only relevant relation is (9). The function
$-f(q)$ decreases from $+\infty$ at $2M^{2}$ to $-1$ at $+\infty$. This means
that $p(q)$ decreases from $+\infty$ at $q=2M^{2}$ until
$q=\frac{m^{4}}{8(m-M)^{2}}$ where $p=0$. This is illustrated on figure 1. The
plot for $p<0$ is just a reflection by the $q$ axis.

Figure 1: Dust shell’s trajectory with energy $H(p,q)=M<\frac{m}{2}$.

For $M=\frac{m}{2}$ the function $f(q)$ takes only values smaller than $1$ and
therefore there is no solution.

If $M>\frac{m}{2}$, we have $\sqrt{q}<\frac{m^{2}}{2\sqrt{2}M}$ so the
function $f(q)$ is fully determined by the relation (8). The case
$M>\frac{m}{2}$ splits further into two other cases:

  1. 1.

$\frac{m}{2}<M<m$: In this case the function $f(q)$ decreases from $+\infty$
at $q=2M^{2}$, reaches a minimal value smaller than $1$ and grows
asymptotically to $1$ as $q$ goes to $+\infty$. This means that $p(q)$
decreases from $+\infty$ at $q=2M^{2}$ until $q=\frac{m^{4}}{8(m-M)^{2}}$
where $p=0$. This is illustrated on figure 2a. The plot for $p<0$ is just a
reflection by the $q$ axis.

  2. 2.

$M\geq m$: In this case the function $f(q)$ decreases from $+\infty$ at
$q=2M^{2}$ to $1$ at $+\infty$. This means that that $p(q)$ decreases from
$+\infty$ at $q=2M^{2}$ to $0$ at $q\rightarrow+\infty$. This is illustrated
on figure 2b.

(a) $\frac{m}{2}<M<m$

(b) $M\geq m$

Figure 2: Dust shells’ trajectories with energy $H(p,q)=M>\frac{m}{2}$.

Let us examine the issue of the positive definiteness of $H(p,q)$. For $p=0$,
Eq. (4) reads $H(0,q)=m-m^{2}/\sqrt{8q}\leavevmode\nobreak\ $ so that we have

$H(0,q)<0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\mbox{for}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
q<m^{2}/8.$ (10)

Since wee need to have $H(p,q)>0$, the physical phase space $\Lambda$ should
be a subspace of the entire phase space $\Pi$. Some regions of $\Pi$, for
instance defined by (10), may lead to the breaking of this condition.

##  III Quantum level

The affine coherent states quantization applied in this article is based on
the formalism presented in our papers [4] and [5].

The physical phase space of our gravitational system

$\Pi:=\\{(p,q)\;|\;p\in\mathbb{R},q\in\mathbb{R}_{+}\\}\,,$ (11)

can be identified with the affine group
$\mathrm{G}\equiv\mathrm{Aff}(\mathbb{R})$, by defining the multiplication law
as follows

$(p^{\prime},q^{\prime})\cdot(p,q):=(q^{\prime}p+p^{\prime},q^{\prime}q)\,,$
(12)

with the unity $(0,1)$ and the inverse

$(p^{\prime},q^{\prime})^{-1}=(-\frac{p^{\prime}}{q^{\prime}},\frac{1}{q^{\prime}})\,.$
(13)

###  III.1 Affine coherent states

The affine group has two, nontrivial, inequivalent irreducible unitary
representations [6, 7, 8]. Both are realized in the Hilbert space
$\mathcal{H}=L^{2}(\mathbb{R}_{+},d\nu(x))$, where $d\nu(x):=dx/x$ is an
invariant measure on the multiplicative group $(\mathbb{R}_{+},\cdot)$. In
what follows we choose the one defined by the following action

$U(p,q)\psi(x)=e^{ipx}\psi(qx)\,,$ (14)

where111We use Dirac’s notation whenever we wish to deal with abstract vector,
instead of functional representation of the vector. $\psi(x)=\langle
x|\psi\rangle$ and $|\psi\rangle\in L^{2}(\mathbb{R}_{+},d\nu(x))$.

We define the integrals over the affine group
$\mathrm{G}=\mathrm{Aff}(\mathbb{R})$ as follows

$\int_{\mathrm{G}}d\mu(p,q)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}dp\int_{0}^{\infty}\frac{dq}{q^{2}}\,.$
(15)

Fixing the normalized vector $|\Phi\rangle\in L^{2}(\mathbb{R}_{+},d\nu(x))$,
called the fiducial vector, we can define a continuous family of affine
coherent states $|p,q\rangle\in L^{2}(\mathbb{R}_{+},d\nu(x))$ as follows

$|p,q\rangle=U(p,q)|\Phi\rangle.$ (16)

The irreducibility of the representation, used to define the coherent states
(16), enables making use of Schur’s lemma [9], which leads to the resolution
of the unity in $L^{2}(\mathbb{R}_{+},d\nu(x))$

$\frac{1}{A_{\Phi}}\int_{\mathrm{G}}d\mu(p,q)|p,q\rangle\langle
p,q|=\mathbb{I}\;,$ (17)

where the constant $A_{\Phi}$ can be determined by using any arbitrary,
normalized vector $|f\rangle\in L^{2}(\mathbb{R}_{+},d\nu(x))$ as follows

$A_{\Phi}=\int_{\mathrm{G}}d\mu(p,q)\,\langle f|p,q\rangle\langle
p,q|f\rangle=\int_{0}^{\infty}\frac{dq}{q^{2}}|\Phi(q)|^{2}\,.$ (18)

###  III.2 Quantum observables

Using the resolution of the identity (17), we define the quantization of a
classical observable $f$ as follows [11]

$\mathcal{F}\ni
f\longrightarrow\hat{f}:=\frac{1}{A_{\Phi}}\int_{\mathrm{G}}d\mu(p,q)|p,q\rangle
f(p,q)\langle p,q|\in\mathcal{A}\,,$ (19)

where $\mathcal{F}$ is a vector space of real continuous functions on a phase
space, and $\mathcal{A}$ is a vector space of operators (quantum observables)
acting in the Hilbert space $L^{2}(\mathbb{R}_{+},d\nu(x))$. It is clear that
(19) defines a linear mapping and the observable $\hat{f}$ is a symmetric
operator. Self-adjointness of $\hat{f}$ is an open problem as symmetricity
does not assure self-adjointness so that further examination is required [12].

In Appendix A we define an orthonormal basis (for any fixed value of the
parameter $\alpha>-1$) of the unitary irreducible representation of considered
affine group. This basis can be used in concrete calculations. To make these
calculations feasible, we use the technics of generating functions for
generalized Laguerre polynomials $L_{n}^{(\alpha)}$ (see, e.g. [13]). For this
purpose, it is convenient to restrict the upper label $\alpha$ of the
functions $L_{n}^{(\alpha)}(x)$ to any fixed integer. In this case, the
generating function for the Laguerre polynomials reads:

$\frac{\exp(-\frac{xz}{1-z})}{(1-z)^{\alpha+1}}=\sum_{n=0}^{\infty}L_{n}^{(\alpha)}(x)z^{n}\quad\mbox{if
}|z|<1\ .$ (20)

For calculation of the matrix elements of the operator (19), in the basis
defined in App. A, one needs to calculate the overlaps between the coherent
state vectors and the vectors of this orthonormal basis. To perform these
calculations explicitly, we choose the fiducial vector for our coherent states
as follows

$\Phi^{(\alpha)}(x):=e^{(\alpha)}_{0}(x)=\sqrt{\frac{1}{\alpha!}}x^{\frac{1+\alpha}{2}}e^{-\frac{x}{2}}\,.$
(21)

In this case one gets

$\langle
e^{(\alpha)}_{n}(x)|pq\rangle=\sqrt{\frac{n!}{(n+\alpha)!\alpha!}}q^{\frac{1+\alpha}{2}}\int_{0}^{\infty}\frac{dx}{x}e^{ipx}x^{1+\alpha}e^{-\frac{q+1}{2}x}L_{n}^{(\alpha)}(x)\,.$
(22)

Let us define the function

$\displaystyle\mathrm{epq}(\alpha,p,q;z):=\sum_{n}\langle
e^{(\alpha)}_{n}(x)|pq\rangle\left[\sqrt{\frac{n!}{(n+\alpha)!\alpha!}}\right]^{-\frac{1}{2}}z^{n}$
$\displaystyle=\left(\frac{\sqrt{q}}{1-z}\right)^{1+\alpha}\int_{0}^{\infty}dx\,x^{\alpha}\exp\left(-\left[\frac{q+1}{2}+\frac{z}{1-z}-ip\right]x\right)\,.$
(23)

The required overlaps turn out to be defined by the derivatives of these
functions calculated at $z=0$ as follows

$\langle
e^{(\alpha)}_{n}(x)|pq\rangle=\frac{1}{\sqrt{\alpha!n!(n+\alpha)!}}\left[\frac{d^{n}}{dz^{n}}\mathrm{epq}(\alpha,p,q;z)\right]_{|z=0}\,.$
(24)

Using exactly the same method, the matrix elements of the operator (19) are
found to be

$\displaystyle\langle
e^{(\alpha)}_{n}|\hat{f}|e^{(\alpha)}_{n^{\prime}}\rangle=\frac{1}{A_{\Phi}}\frac{1}{\sqrt{\alpha!n!(n+\alpha)!\alpha!{n^{\prime}}!(n^{\prime}+\alpha)!}}$
$\displaystyle\cdot\left[\frac{\partial^{n}}{\partial
z_{1}^{n}}\frac{\partial^{n^{\prime}}}{\partial
z_{2}^{n^{\prime}}}\int_{\mathrm{Aff}(\mathbb{R})}d\mu_{L}(p,q)\mathrm{epq}(\alpha,p,q;z_{1})f(p,q)\mathrm{epq}(\alpha,p,q;z_{2})^{\star}\right]_{|z_{1}=z_{2}=0}\,.$
(25)

The simplest basis, satisfying required conditions, is obtained by taking
$\alpha=1$. In this case, the function $\mathrm{epq}(1,p,q;z)$ reads:

$\mathrm{epq}(1,p,q;z)=\frac{q}{(1-z)^{2}}\left(\frac{q+1}{2}+\frac{z}{1-z}-ip\right)^{-2}\,.$
(26)

###  III.3 Quantum dynamics

The mapping (19) applied to the classical Hamiltonian (4) reads

$\hat{H}_{\rm unbounded}=\frac{1}{A_{\Phi}}\int_{G}d\mu(p,q)|p,q\rangle
H(p,q)\langle p,q|\,.$ (27)

However, our classical analysis applies only to the region of phase space for
which $H(p,q)>0$.

An important problem is introducing constraints into the integral quantization
approach. This quantization is based on deformation of quantum measure
represented by a set of positive self-adjoint operators determining the
operator valued measure (POV). They are considered as generalization of more
standard self-adjoint quantum observables [14].

In our case the set of operators

$\hat{M}(Q):=\frac{1}{A_{\phi}}\int_{G}d\mu(p,q)|p,q\rangle\chi_{Q}(g)\langle
p,q|\,,$ (28)

where $\chi_{Q}(g)=1$ if $g\in Q$ and $0$ otherwise, and where
$Q\subset\mathrm{G}=\mathrm{Aff}(\mathbb{R})$, describe the localization of
the system in the subspace $Q$ of the phase space $\mathrm{G}$.

The normalization condition

$\hat{M}(G)=\frac{1}{A_{\phi}}\int_{G}d\mu(p,q)|p,q\rangle\langle
p,q|=\hat{1\kern-4.75pt1}\,$ (29)

is required to get the so called minimal probabilistic interpretation of
quantum mechanics [14]:

$\mathrm{Prob}(Q,\Psi)=\langle\Psi|\hat{M}(\mathrm{Q})|\kern
1.29167pt\Psi\kern 1.29167pt\rangle\,,$ (30)

which describes the probability of finding our system in $Q$, under condition
that this system is in the state $\Psi$.

On the other hand, to construct the condition (29) one needs to integrate over
the whole group manifold and the representation (16) has to be irreducible.
This excludes the possibility of using a smaller region
$\Omega\subset\mathrm{G}$ of the phase space to fulfill the classical
constraint, where $\Omega=\\{(p,q):H(p,q)>0\\}$.

To have a consistent quantization method, the only possibility of quantizing
any observable restricted to a smaller region of the phase space is to
quantize it over the whole phase space, represented in this approach by the
group $\mathrm{G}$.

Following these requirements, to keep physical interpretation of the classical
Hamiltonian, we quantize it as a function restricted to the required region of
the phase space by considering:

$\hat{H}=\frac{1}{A_{\Phi}}\int_{G}d\mu(p,q)|p,q\rangle\theta(H(p,q))H(p,q)\langle
p,q|\,,$ (31)

where $\theta$ is the Heaviside theta function.

An important feature of the operator $\hat{H}$ is that it acts in a nontrivial
way on the whole phase space. Let $(p^{\prime},q^{\prime})\not\in\Omega$, then

$\hat{H}|p^{\prime},q^{\prime}\rangle=\frac{1}{A_{\Phi}}\int_{\Omega}d\mu(p,q)|p,q\rangle\theta(H(p,q))H(p,q)\langle
p,q|p^{\prime},q^{\prime}\rangle$ (32)

is usually a non-zero vector. This behavior is due to non-orthogonality of the
states $|p,q\rangle$, i.e., the quantum phase space regions $\Omega$ and
$\mathrm{G}\setminus\Omega$ are not independent.

Another very important problem related is the quantization of the elementary
observables $(p,q)$. These observables represent position of the system in the
quantum full phase space $\mathrm{G}$. The corresponding quantum operators
$(\hat{p},\hat{q})$ should satisfy the following consistency conditions [15]:

$\langle p,q|\hat{p}|p,q\rangle=p\text{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ }\langle
p,q|\hat{q}|p,q\rangle=q\,,$ (33)

where

$\hat{p}=\frac{1}{A_{\Phi}}\int_{G}d\mu(p,q)|p,q\rangle p\langle
p,q|\text{\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
and\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ }\hat{q}=\frac{1}{A_{\Phi}}\int_{G}d\mu(p,q)|p,q\rangle
q\langle p,q|\,.$ (34)

Satisfying (33) is possible only when the operators are defined on the entitre
phase space $\mathrm{G}$. The consistency conditions support the physical
interpretation of the POV measure (28).

The above analysis allow us to be consistent and to quantize the restricted
form of the Hamiltonian $H$ and unrestricted form of the elementary
observables $(p,q)$.

We will restrict to this region of a phase space by considering:

$\hat{H}=\frac{1}{A_{\Phi}}\int_{G}d\mu(p,q)|p,q\rangle\theta(H(p,q))H(p,q)\langle
p,q|\,,$ (35)

where $\theta$ is the Heaviside theta function. Let us notice that the
operator $\hat{H}$ is positive definite. Indeed, for any state $|\kern
1.29167pt\Psi\kern 1.29167pt\rangle$ we have:

$\langle\kern 1.29167pt\Psi\kern 1.29167pt|\kern 1.29167pt\hat{H}\Psi\kern
1.29167pt\rangle=\frac{1}{A_{\Phi}}\int_{G}d\mu(p,q)|\Psi(p,q)|^{2}\theta(H(p,q))H(p,q)\,>0.$
(36)

The quantum evolution of our gravitational system is defined by the
Schrödinger equation:

$i\frac{\partial}{\partial s}|\Psi(s)\rangle=\hat{H}|\Psi(s)\rangle\;,$ (37)

where $|\Psi\rangle\in L^{2}(\mathbb{R}_{+},d\nu(x))$, and where $s$ is an
evolution parameter of the quantum level.

In general, the parameters $t$ of the classical level and $s$ are different.
To get the consistency between the classical and quantum levels we postulate
that $t=s$, which defines the time variable at both levels. This way we
support the interpretation that Hamiltonian is the generator of classical and
corresponding quantum dynamics.

##  IV Matrix elements of quantum observables

In this section we calculate the matrix elements for the operators
$\hat{H},\hat{q}$, and $\hat{p}$. The matrix elements of $\hat{H}$ are
calculated numerically using (60). We calculate the matrix elements of
$\hat{q}$ and $\hat{p}$ analytically. The computations are based on a new
formula for the basis elements $e^{(1)}_{m}$ that we derive in the first
subsection. The results of this section are used in Sec. V to find an
evolution of the quantum observables.

###  IV.1 New expression for the basis elements

In what follows, we use the basis with $\alpha=1$. Since
$\frac{z}{1-z}=\frac{1}{1-z}-1$ and
$\frac{d}{dz}\frac{1}{1-z}=\frac{1}{(1-z)^{2}}$, we can simplify the
expression for the basis element:

$\langle
e^{(1)}_{n}|pq\rangle=-\frac{q\sqrt{n+1}}{(n+1)!}\left[\frac{d^{n+1}}{dz^{n+1}}\int_{0}^{\infty}dx\,\exp\left(-\left[\frac{q-1}{2}+\frac{1}{1-z}-ip\right]x\right)\right]_{|z=0}$
(38)

The technical problem is now to calculate $\frac{d^{n}}{dz^{n}}f(z)$, where
$f(z)$ is a composite function:

$f(z)=F(G(z))=\exp\left(-G(z)x\right),{\rm\ where\ }G(z)=\frac{1}{1-z}.$ (39)

We will use the Faa di Bruno’s formula [16, 17]:

$\frac{d^{n}}{dz^{n}}F(G(z))=\sum_{k=1}^{n}F^{(k)}(G(z))\cdot
B_{n,k}(G^{\prime}(z),G^{\prime\prime}(z),\ldots,G^{(n-k+1)}(z)),$ (40)

where $B_{n,k}(x_{1},\ldots,x_{n-k+1})$ are Bell polynomials [18, 17]. In the
formula above we denote by $F^{(k)}$ the k-th derivative of $F$. The Bell
polynomial $B_{n,k}(z_{1},z_{2},\ldots,z_{n-k+1})$ is given by

$B_{n,k}(z_{1},z_{2},\ldots,z_{n-k+1})=\sum\frac{n!}{j_{1}!j_{2}!\ldots
j_{n-k+1}!}\left(\frac{z_{1}}{1!}\right)^{j_{1}}\left(\frac{z_{2}}{2!}\right)^{j_{2}}\ldots\left(\frac{z_{n-k+1}}{(n-k+1)!}\right)^{j_{n-k+1}},$
(41)

where the sum is over all sequences $j_{1},\ldots,j_{n-k+1}\in\mathbb{N}$ such
that

$\displaystyle j_{1}+j_{2}+\ldots+j_{n-k+1}=k,$ (42) $\displaystyle
j_{1}+2\,j_{2}+3\,j_{3}+\ldots+(n-k+1)j_{n-k+1}=n.$ (43)

We are interested only in derivatives at $z=0$, therefore in equation (40) we
should put $z=0$. Let us notice that

$G(0)=1,\quad G^{m}(0)=m!\,.$ (44)

Therefore

$\frac{d^{n}}{dz^{n}}F(G(0))=\sum_{k=1}^{n}(-1)^{k}x^{k}e^{-x}\,B_{n,k}(1!,2!,\ldots,(n-k+1)!).$
(45)

The coefficient $B_{n,k}(1!,2!,\ldots,(n-k+1)!)$ can be expressed in terms of
the Lah number [17, 19]:

$B_{n,k}(1!,2!,\ldots,(n-k+1)!)=|L_{n,k}|,$ (46)

where

$\displaystyle L_{n,k}=(-1)^{n}\binom{n-1}{k-1}\frac{n!}{k!},\quad n\geq k\geq
1,$ (47) $\displaystyle L_{0,0}=1,\quad L_{n,0}=0,\,n\geq 1.$ (48)

We extend the sum in (45) to $k=0$ and obtain a formula valid for $n=0$:

$\frac{d^{n}}{dz^{n}}F(G(0))=\sum_{k=0}^{n}(-1)^{k}|L_{n,k}|\,x^{k}e^{-x}.$
(49)

In order to calculate $\langle e^{(1)}_{n}|pq\rangle$, we need to evaluate the
integral

$\int_{0}^{\infty}dx\
x^{k}\exp\left(-\left(\frac{q+1}{2}-\texttt{i}p\right)x\right)=\frac{k!}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{k+1}}.$
(50)

This integral directly follows from the Euler’s integral formula (see for
example equation 6.1.1 from [20]:

$\Gamma(n)=z^{n}\int_{0}^{\infty}dx\,x^{n-1}e^{-zx},\leavevmode\nobreak\
\quad{\rm for\ }\leavevmode\nobreak\ \Re(n)>0,\leavevmode\nobreak\
\Re{(z)}>0.$ (51)

In the formula above $\Gamma(z)$ is the Gamma function. Inserting this
evaluation into (38) gives:

$\langle
e^{(1)}_{n}|pq\rangle=\sqrt{n+1}\sum_{k=0}^{n+1}\frac{(-1)^{k+1}k!}{(n+1)!}|L_{n+1,k}|\frac{q}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{k+1}}\,.$
(52)

We will define coefficients $E_{n,k}$:

$E_{n,k}=\frac{(-1)^{k+1}k!}{n!}|L_{n,k}|.$ (53)

Due to cancellation of terms, the coefficients have the following explicit
form:

$\displaystyle E_{n,k}=(-1)^{k-1}\binom{n-1}{k-1},\quad n\geq k\geq 1,$ (54)
$\displaystyle E_{0,0}=\frac{(-1)^{k-1}k!}{n!},\quad E_{n,0}=0,\,n\geq 1.$
(55)

Let us notice, that in our formulas we need only the coefficients with $n\geq
1$. Since, $E_{n,0}=0$ for $n\geq 1$, we can assume that $k\geq 1$. The
formula (52) takes now a more compact form:

$\langle
e^{(1)}_{n}|pq\rangle=\sqrt{n+1}\sum_{k=1}^{n+1}(-1)^{k-1}\binom{n}{k-1}\frac{q}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{k+1}}\,.$
(56)

We will change the summation variable into $k^{\prime}=k-1$ and obtain:

$\langle
e^{(1)}_{n}|pq\rangle=\frac{q\,\sqrt{n+1}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{2}}\,\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{k}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{k}}.$
(57)

The last sum is the binomial expansion. Therefore

$\langle
e^{(1)}_{n}|pq\rangle=\frac{q\,\sqrt{n+1}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{2}}\,\left(1-\frac{1}{\frac{q+1}{2}-\texttt{i}p}\right)^{n}=q\,\sqrt{n+1}\,\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{n+2}}.$
(58)

###  IV.2 The matrix elements of the Hamiltonian operator

The new expression for the basis elements that we found in the previous
section allows us to calculate the matrix elements of the Hamiltonian:

$H_{nm}:=\langle
e^{(1)}_{n}|\hat{H}e^{(1)}_{m}\rangle=\frac{1}{A_{\Phi}}\int_{G}d\mu(p,q)\langle
e^{(1)}_{n}|pq\rangle\theta(H(q,p))H(p,q)\langle e^{(1)}_{m}|pq\rangle^{*}.$
(59)

After applying the formula (58) to the equation above, we obtain:

$H_{nm}:=\frac{\sqrt{(n+1)(m+1)}}{2\pi}\,\int_{0}^{\infty}dq\int_{-\infty}^{+\infty}dp\
\theta(H(q,p))H(q,p)\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{n+2}\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}.$
(60)

Since $\hat{H}$ is hermitian, it is enough to calculate the lower triangular
part of the matrix, i.e. $H_{nm}$ for $n\geq m$. In this case:

$H_{nm}=\frac{\sqrt{(n+1)(m+1)}}{2\pi}\int_{0}^{\infty}dq\int_{-\infty}^{+\infty}dp\
\theta(H(q,p))H(q,p)\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}\left(\frac{q+1}{2}+\texttt{i}p\right)^{n-m}}{\left(\left(\frac{q+1}{2}\right)^{2}+p^{2}\right)^{n+2}}.$
(61)

In order to calculate the real and the imaginary part of the matrix, we write
the Hamiltonian in the form

$H_{nm}=\frac{\sqrt{(n+1)(m+1)}}{2\pi}\int_{0}^{\infty}dq\int_{-\infty}^{+\infty}dp\
\theta(H(q,p))H(q,p)\frac{\left(\left(\frac{q-1}{2}\right)^{2}+p^{2}\right)^{m}\left(\frac{q^{2}-1}{4}+p^{2}-\texttt{i}p\right)^{n-m}}{\left(\left(\frac{q+1}{2}\right)^{2}+p^{2}\right)^{n+2}}.$
(62)

After expanding the expression
$\left(\frac{q^{2}-1}{4}+p^{2}-\texttt{i}p\right)^{n-m}$ in powers of i, we
notice that the real part of the integrant is an even function in $p$ and the
imaginary part is an odd function in $p$. As a result, the imaginary part of
$H_{nm}$ vanishes:

$\Im H_{nm}=0.$ (63)

We conclude that $H_{nm}$ is a real symmetric matrix. We will calculate the
matrix elements by evaluating the integrals over $q$ and $p$ numerically.

The diagonal elements $H_{nn}$ are positive since the integrand of (62) is a
positive function of the variables $q$ and $p$.

###  IV.3 The matrix elements of the momentum operator

The matrix elements of the momentum operator:

$p_{nm}:=\frac{\sqrt{(n+1)(m+1)}}{2\pi}\,\int_{0}^{\infty}dq\int_{-\infty}^{+\infty}dp\
\frac{p\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{n+2}\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}.$
(64)

The integral over $p$ can be done using a contour method. We take the sunset
contour in the lower half plane and obtain

$p_{nm}:=-\texttt{i}\sqrt{(n+1)(m+1)}\,\int_{0}^{\infty}dq\ {\rm
Res}|_{p=-\texttt{i}\frac{q+1}{2}}\left(\frac{p\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{n+2}\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}\right).$
(65)

We will calculate now the residuum:

${\rm
Res}|_{p=-\texttt{i}\frac{q+1}{2}}\left(\frac{p\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{n+2}\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}\right)=\\\
=\texttt{i}^{n+2}\frac{1}{(n+1)!}\frac{d^{n+1}}{dp^{n+1}}\left({\frac{p\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}}\right)=\\\
=\frac{\texttt{i}^{n+2}}{(n+1)!}p\frac{d^{n+1}}{dp^{n+1}}\left({\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}}\right)|_{p=-\texttt{i}\frac{q+1}{2}}+\\\
+\frac{\texttt{i}^{n+2}}{n!}\frac{d^{n}}{dp^{n}}\left({\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}}\right)|_{p=-\texttt{i}\frac{q+1}{2}}.$
(66)

Let us notice that the two terms in the last line are proportional to
expressions of the form:

$I(l,n,m,q)=\frac{1}{l!}\frac{d^{l}}{dp^{l}}\left({\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}}\right)|_{p=-\texttt{i}\frac{q+1}{2}}$
(67)

with $l=n$ or $l=n+1$. We will expand the derivative using a Leibnitz rule
generalized to three factors:

$I(l,n,m,q)=\frac{1}{l!}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=l\end{subarray}}\binom{l}{k_{1},k_{2},k_{3}}\frac{d^{k_{1}}}{dp^{k_{1}}}\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\cdot\\\
\cdot\frac{d^{k_{2}}}{dp^{k_{2}}}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}\frac{d^{k_{3}}}{dp^{k_{3}}}\left(\frac{q+1}{2}+\texttt{i}p\right)^{-m-2}|_{p=-\texttt{i}\frac{q+1}{2}},$
(68)

where $\binom{l}{k_{1},k_{2},k_{3}}=\frac{l!}{k_{1}!k_{2}!k_{3}!}$. The
derivatives of power functions can be calculated and give
($j,k\in\mathbb{Z},j\geq 0,\ k\geq 0$):

$\displaystyle\frac{d^{k}}{dx^{k}}x^{j}=\begin{cases}\frac{j!}{(j-k)!}\cdot
x^{j-k},&{\rm for}\ j\geq k,\\\ 0,&{\rm for}\ j<k,\end{cases}$ (69)
$\displaystyle\frac{d^{k}}{dx^{k}}x^{-j}=(-1)^{k}\frac{(j+k-1)!}{(j-1)!}\cdot
x^{-j-k}.$ (70)

Inserting this result into (68) gives:

$I(l,n,m,q)=\frac{1}{l!}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=l\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}\binom{l}{k_{1},k_{2},k_{3}}\frac{n!}{(n-k_{1})!}(-\texttt{i})^{k_{1}}\left(\frac{q-1}{2}-\texttt{i}p\right)^{n-k_{1}}\cdot\\\
\cdot(\texttt{i})^{k_{2}}\frac{m!}{(m-k_{2})!}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m-k_{2}}(-\texttt{i})^{k_{3}}\frac{(m+1+k_{3})!}{(m+1)!}\left(\frac{q+1}{2}+\texttt{i}p\right)^{-m-2-k_{3}}|_{p=-\texttt{i}\frac{q+1}{2}}.$
(71)

After evaluating at $p=-\texttt{i}\frac{q+1}{2}$ we get:

$I(l,n,m,q)=(-1)^{n}(\texttt{i})^{l}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=l\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\binom{n}{k_{1}}\binom{m}{k_{2}}\binom{m+1+k_{3}}{k_{3}}\frac{q^{m-k_{2}}}{\left(q+1\right)^{m+k_{3}+2}}.$
(72)

In order to find the expression for the matrix elements of the momentum
operator, we need to evaluate 2 integrals:

$I_{1}(n,m)=\int_{0}^{\infty}dq\ (q+1)I(n+1,n,m,q),\quad
I_{2}(n,m)=\int_{0}^{\infty}dq\ I(n,n,m,q).$ (73)

In order to perform the integrals, we notice that they can be expressed in
terms of the beta function

$B(x,y)=\int_{0}^{\infty}dt\frac{t^{x-1}}{(1+t)^{x+y}},\quad{\rm for\
}\leavevmode\nobreak\ \Re{(x)}>0,\leavevmode\nobreak\ \Re{y}>0.$ (74)

The first integral is:

$I_{1}(n,m)=(-1)^{n}(\texttt{i})^{n+1}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\binom{n}{k_{1}}\binom{m}{k_{2}}\binom{m+1+k_{3}}{k_{3}}B(m-k_{2}+1,k_{2}+k_{3}).$
(75)

If $x,y\in\mathbb{N}-\\{0\\}$, the beta function takes the form:

$B(x,y)=\frac{(x-1)!(y-1)!}{(x+y-1)!}.$ (76)

Inserting this property into equation (75) we get:

$I_{1}(n,m)=(-1)^{n}\texttt{i}^{n+1}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\binom{n}{k_{1}}\binom{m}{k_{2}}\binom{m+1+k_{3}}{k_{3}}\frac{(m-k_{2})!(k_{2}+k_{3}-1)!}{(m+k_{3})!}=\\\
=\frac{(-1)^{n}\texttt{i}^{n+1}}{(n+1)(m+1)}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\frac{(n+1)!}{k_{1}!k_{2}!k_{3}!}\left(m+1+k_{3}\right)=\\\
=\frac{(-1)^{n}\texttt{i}^{n+1}}{(n+1)(m+1)}\left(-(m+1)+\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n+1,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\frac{(n+1)!}{k_{1}!k_{2}!k_{3}!}\left(m+1+k_{3}\right)\right)\,.$
(77)

The second integral is:

$I_{2}(n,m)=(-1)^{n}(\texttt{i})^{n}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\binom{n}{k_{1}}\binom{m}{k_{2}}\binom{m+1+k_{3}}{k_{3}}\frac{(m-k_{2})!(k_{2}+k_{3})!}{(m+k_{3}+1)!}=\\\
=\frac{(-1)^{n}(\texttt{i})^{n}}{m+1}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\frac{n!}{k_{1}!k_{2}!k_{3}!}\,.$ (78)

The matrix elements of the momentum operator can be expressed in terms of the
functions $I_{1}(n,m)$ and $I_{2}(n,m)$:

$p_{nm}:=\sqrt{(n+1)(m+1)}\texttt{i}^{n+2}\left(-\frac{1}{2}I_{1}(n,m)-\texttt{i}I_{2}(n,m)\right)$
(79)

####  IV.3.1 Entries below and on the diagonal

We will look for a compact formula for the matrix elements. To this end we
will find expressions for

$S_{n,m}=\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\ k_{1}+k_{2}+k_{3}=n\\\
k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\frac{n!}{k_{1}!k_{2}!k_{3}!},\quad
Z_{n,m}=\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\ k_{1}+k_{2}+k_{3}=n\\\
k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}k_{3}\frac{n!}{k_{1}!k_{2}!k_{3}!}.$ (80)

We will find $S_{n,m}$ and $Z_{n,m}$ recursively. In this derivation we will
rely on the following identity holding for any polynomial $P$ of degree
smaller than $n$:

$\sum_{k=0}^{n}(-1)^{k}\binom{n}{k}P(k)=0.$ (81)

Let us consider first $S_{n,m}$. From the definition, it is clear that

$S_{n,m+1}=S_{n,m}+\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n\\\ k_{1}\leq
n,k_{2}=m+1\end{subarray}}(-1)^{k_{3}}\frac{n!}{k_{1}!k_{2}!k_{3}!}=\\\
=S_{n,m}+\frac{n!}{(n-m-1)!(m+1)!}\sum_{\begin{subarray}{c}k_{1},k_{3}\\\
k_{1}+k_{3}=n-m-1\end{subarray}}(-1)^{k_{3}}\frac{(n-m-1)!}{k_{1}!k_{3}!}=\\\
=\begin{cases}S_{n,m}+1,&\textrm{ if }m=n-1,\\\ S_{n,m},&\textrm{ if
}m<n-1.\\\ \end{cases}$ (82)

In order to find the desired identity, we notice that

$S_{n,0}=\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\ k_{1}+k_{2}+k_{3}=n\\\
k_{1}\leq
n,k_{2}=0\end{subarray}}(-1)^{k_{3}}\frac{n!}{k_{1}!k_{2}!k_{3}!}=\sum_{\begin{subarray}{c}k_{1},k_{3}\\\
k_{1}+k_{3}=n\end{subarray}}(-1)^{k_{3}}\frac{n!}{k_{1}!k_{3}!}=\delta_{n,0}.$
(83)

Combining the results we get that :

$S_{n,m}=\delta_{n,m}\textrm{ for }m\leq n.$ (84)

The second sum will be calculated in completely similar manner. We have the
following recurrence:

$Z_{n,m+1}=Z_{n,m}+\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n\\\ k_{1}\leq
n,k_{2}=m+1\end{subarray}}(-1)^{k_{3}}k_{3}\,\frac{n!}{k_{1}!k_{2}!k_{3}!}=\\\
=Z_{n,m}+\binom{n}{m+1}\sum_{\begin{subarray}{c}k_{1},k_{3}\\\
k_{1}+k_{3}=n-m-1\end{subarray}}(-1)^{k_{3}}k_{3}\binom{n-m-1}{k_{3}}=\\\
=\begin{cases}Z_{n,m},&\textrm{ if }m=n-1,\\\ Z_{n,m}-n,&\textrm{ if
}m=n-2.\\\ Z_{n,m},&\textrm{ if }m\leq n-3.\end{cases}$ (85)

We can express the formula above using Kronecker delta:

$Z_{n,m+1}=Z_{n,m}-n\delta_{n,m+2}.$ (86)

The recurrence starts at $Z_{n,0}$:

$Z_{n,0}=\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\ k_{1}+k_{2}+k_{3}=n\\\
k_{1}\leq
n,k_{2}=0\end{subarray}}(-1)^{k_{3}}k_{3}\,\frac{n!}{k_{1}!k_{2}!k_{3}!}=\sum_{\begin{subarray}{c}k_{1},k_{3}\\\
k_{1}+k_{3}=n\end{subarray}}(-1)^{k_{3}}k_{3}\,\binom{n}{k_{3}}=\begin{cases}1,&\textrm{
if }n=0,\\\ -1,&\textrm{ if }n=1.\\\ 0,&\textrm{ if }n\geq 2.\end{cases}$ (87)

In our study, we will only need to consider the following two cases.

  1. 1.

$n=1$. In this case $Z_{1,0}=-1$ and $Z_{1,1}=-1$.

  2. 2.

$n\geq 2$. In this case $Z_{n,0}=0$ and

$Z_{n,m+1}=Z_{n,m}-n\delta_{n,m+2}.$ (88)

As a result,

$Z_{n,m}=-n\delta_{n,m+1}-n\delta_{n,m}.$ (89)

Let us notice that the two cases can be written with one formula. For $n\geq
1$:

$Z_{n,m}=-n\delta_{n,m+1}-n\delta_{n,m}.$ (90)

Let us now go back to the expression for $I_{1}(n,m)$ and $I_{2}(n,m)$.

$I_{1}(n,m)=\frac{(-1)^{n}\texttt{i}^{n+1}}{(n+1)(m+1)}\left((m+1)(S_{n+1,m}-1)+Z_{n+1,m}\right)=\\\
=\frac{(-1)^{n}\texttt{i}^{n+1}}{(n+1)(m+1)}\left((m+1)(\delta_{n+1,m}-1)-(n+1)\delta_{n+1,m+1}-(n+1\delta_{n+1,m}\right).$
(91)

The second object is:

$I_{2}(n,m)=\frac{(-1)^{n}(\texttt{i})^{n}}{m+1}S_{n,m}=\frac{(-1)^{n}(\texttt{i})^{n}}{m+1}\delta_{n,m}.$
(92)

We will combine the results to obtain the expression for the matrix elements
of the momentum operator:

$p_{nm}=\frac{\texttt{i}}{\sqrt{(n+1)(m+1)}}\left(\frac{m+1}{2}(\delta_{n+1,m}-1)-\frac{n+1}{2}\delta_{n+1,m+1}-\frac{n+1}{2}\delta_{n+1,m}+(n+1)\delta_{n,m}\right).$
(93)
$p_{nm}=\frac{\texttt{i}}{\sqrt{(n+1)(m+1)}}\left(-\frac{m+1}{2}+\frac{n+1}{2}\delta_{n,m}+\frac{m-n}{2}\delta_{n+1,m}\right).$
(94)

Let us notice that we make calculation for $m\leq n$ only. Therefore

$p_{nm}=\frac{\texttt{i}}{\sqrt{(n+1)(m+1)}}\left(-\frac{m+1}{2}+\frac{n+1}{2}\delta_{n,m}\right).$
(95)

This means that:

$p_{nn}=0,\quad p_{nm}=-\frac{\texttt{i}}{2}\sqrt{\frac{m+1}{n+1}}\textrm{ if
}m<n.$ (96)

####  IV.3.2 The entries above the diagonal

For a cross-check, we will calculate the entries above the diagonal as well
$m>n$. In this case

$S_{n,m}=\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\ k_{1}+k_{2}+k_{3}=n\\\
k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\frac{n!}{k_{1}!k_{2}!k_{3}!}=\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n\\\ k_{1}\leq n,k_{2}\leq
n\end{subarray}}(-1)^{k_{3}}\frac{n!}{k_{1}!k_{2}!k_{3}!}=(1+1-1)^{n}=1$ (97)

and

$Z_{n,m}=\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\ k_{1}+k_{2}+k_{3}=n\\\
k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}k_{3}\frac{n!}{k_{1}!k_{2}!k_{3}!}=\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n\\\ k_{1}\leq n,k_{2}\leq
n\end{subarray}}(-1)^{k_{3}}k_{3}\frac{n!}{k_{1}!k_{2}!k_{3}!}=\\\
=x_{3}\frac{\partial}{\partial
x_{3}}(x_{1}+x_{2}-x_{3})^{n}|_{x_{1}=x_{2}=x_{3}=1}=-n(1+1-1)^{n}=-n.$ (98)

Inserting the formulas into $I_{1}(n,m),I_{2}(n,m)$ and into $p_{nm}$
afterwards, we get:

$p_{nm}=\frac{\texttt{i}}{2}\sqrt{\frac{n+1}{m+1}}\textrm{ if }m>n.$ (99)

This is consistent with the fact that the matrix $p$ is hermitian.

####  IV.3.3 Summary

The calculations above show that

$p_{nn}=0,\quad p_{nm}=-\frac{\texttt{i}}{2}\sqrt{\frac{m+1}{n+1}}\textrm{ if
}m<n,\quad p_{nm}=\frac{\texttt{i}}{2}\sqrt{\frac{n+1}{m+1}}\textrm{ if }m>n.$
(100)

###  IV.4 The matrix elements of the coordinate operator

The matrix elements of the coordinate operator:

$q_{nm}:=\frac{\sqrt{(n+1)(m+1)}}{2\pi}\,\int_{0}^{\infty}dq\int_{-\infty}^{+\infty}dp\
\frac{q\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{n+2}\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}.$
(101)

The integral over $p$ can be done using a contour method. We take the sunset
contour in the lower half plane and obtain

$q_{nm}:=-\texttt{i}\sqrt{(n+1)(m+1)}\,\int_{0}^{\infty}dq\ q\,{\rm
Res}|_{p=-\texttt{i}\frac{q+1}{2}}\left(\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{n+2}\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}\right).$
(102)

We will calculate now the residuum:

${\rm
Res}|_{p=-\texttt{i}\frac{q+1}{2}}\left(\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}-\texttt{i}p\right)^{n+2}\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}\right)=\\\
=\texttt{i}^{n+2}\frac{1}{(n+1)!}\frac{d^{n+1}}{dp^{n+1}}\left({\frac{\left(\frac{q-1}{2}-\texttt{i}p\right)^{n}\left(\frac{q-1}{2}+\texttt{i}p\right)^{m}}{\left(\frac{q+1}{2}+\texttt{i}p\right)^{m+2}}}\right).$
(103)

The integral that we will need to calculate this time is

$I_{3}(n,m)=\int_{0}^{\infty}dq\,q\,I(n+1,n,m,q).$ (104)

With this integral, we can express the matrix elements of the coordinate
operator in the following form:

$q_{nm}=\texttt{i}^{n+1}\sqrt{(n+1)(m+1)}I_{3}(n,m).$ (105)

As previously, we can express the integral $I_{3}(n,m)$ as a sum.

$I_{3}(n,m)=(-1)^{n}\texttt{i}^{n+1}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\binom{n}{k_{1}}\binom{m}{k_{2}}\binom{m+1+k_{3}}{k_{3}}\int_{0}^{\infty}dq\,\frac{q^{m-k_{2}+1}}{(q+1)^{m+k_{3}+2}}=\\\
=(-1)^{n}\texttt{i}^{n+1}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\binom{n}{k_{1}}\binom{m}{k_{2}}\binom{m+1+k_{3}}{k_{3}}\frac{(m-k_{2}+1)!(k_{2}+k_{3}-1)!}{(m+k_{3}+1)!}=\\\
=-\frac{(-\texttt{i})^{n+1}}{(n+1)(m+1)}\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\frac{(n+1)!}{k_{1}!k_{2}!k_{3}!}(m+1-k_{2})=\\\
=\frac{(-\texttt{i})^{n+1}}{(n+1)(m+1)}\left((m+1)-\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n+1,k_{2}\leq
m\end{subarray}}(-1)^{k_{3}}\frac{(n+1)!}{k_{1}!k_{2}!k_{3}!}(m+1-k_{2})\right).$
(106)

This time we will consider separately the terms on the diagonal and above the
diagonal (and we will not calculate the terms below diagonal.

####  IV.4.1 Diagonal terms

$q_{nn}=\frac{1}{n+1}\left((n+1)-\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n+1,k_{2}\leq
n\end{subarray}}(-1)^{k_{3}}\frac{(n+1)!}{k_{1}!k_{2}!k_{3}!}(n+1-k_{2})\right)=\\\
=\frac{1}{n+1}\left((n+1)-\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n+1,k_{2}\leq
n+1\end{subarray}}(-1)^{k_{3}}\frac{(n+1)!}{k_{1}!k_{2}!k_{3}!}(n+1-k_{2})\right)=\\\
=\frac{1}{n+1}\left((n+1)-(n+1)(1+1-1)^{n+1}+x_{2}\frac{\partial}{\partial
x_{2}}\left(x_{1}+x_{2}-x_{3}\right)^{n+1}|_{x_{1}=x_{2}=x_{3}=1}\right)=1.$
(107)

####  IV.4.2 The terms above the diagonal

We will consider matrix elements $q_{nm}$ for $m>n$. In this case:

$q_{nm}=\frac{1}{\sqrt{(n+1)(m+1)}}\left((m+1)-\sum_{\begin{subarray}{c}k_{1},k_{2},k_{3}\\\
k_{1}+k_{2}+k_{3}=n+1\\\ k_{1}\leq n+1,k_{2}\leq
n+1\end{subarray}}(-1)^{k_{3}}\frac{(n+1)!}{k_{1}!k_{2}!k_{3}!}(m+1-k_{2})\right)=\\\
=\frac{1}{\sqrt{(n+1)(m+1)}}\left((m+1)-(m+1)(1+1-1)^{n+1}+x_{2}\frac{\partial}{\partial
x_{2}}\left(x_{1}+x_{2}-x_{3}\right)^{n+1}|_{x_{1}=x_{2}=x_{3}=1}\right)=\\\
=\sqrt{\frac{n+1}{m+1}}.$ (108)

####  IV.4.3 Summary

Since $q_{nm}$ is hermitian, we have

$q_{nm}=\sqrt{\frac{n+1}{m+1}}\textrm{ if }m\geq n,\quad
q_{nm}=\sqrt{\frac{m+1}{n+1}}\textrm{ if }m<n.$ (109)

##  V Evolution of quantum observables

In what follows, we consider the quantum evolution corresponding to the
classical case with $M<\frac{m}{2}$, where $m$ is the total rest mass of the
shell. Classically, the shell reaches the horizon $q_{H}=2M^{2}$ after an
infinite time (of an observer at spatial infinity), as it is illustrated in
Fig. 1.

Remarkably, our quantum model experiences a different behaviour. In a finite
time, the quantum shell reaches a minimum size and after a bounce, which is
above the horizon, it expands until it reaches a maximum size. More precisely,
in our quantum model we choose a state peaked at energy $M$. One needs to
remember that the energy of the system is conserved during the time evolution
given by the Schrödinger equation. We diagonalize the quantum Hamiltonian
numerically. Due to technical limitations of our numerical procedure we can
study the model for small values of $M$ only. We calculate the quantum
evolution of the expectation values of the operators $\hat{q}$ and $\hat{p}$.
We notice that the expectation values $\langle q\rangle_{t}$ and $\langle
p\rangle_{t}$ are periodic functions of time $t$.

###  V.1 Gaussian state

We choose a state peaked at energy $M$ with standard deviation $\sigma$:

$|\kern 1.29167pt\Psi\kern
1.29167pt\rangle=\frac{1}{N}\sum_{i}e^{-\frac{(E_{i}-M)^{2}}{2\sigma^{2}}}|\kern
1.29167ptE_{i}\kern 1.29167pt\rangle,$ (110)

where $|\kern 1.29167ptE_{i}\kern 1.29167pt\rangle$ is the eigenstate of the
Hamiltonian operator with the eigenvalue $E_{i}$ and $N$ is the normalization
constant:

$N^{2}=\sum_{i}e^{-\frac{(E_{i}-M)^{2}}{\sigma^{2}}}.$ (111)

In practical calculations we limit the sum to the region
$[M-4\sigma\sqrt{2\ln(10)},M+4\sigma\sqrt{2\ln(10)}]$. With this choice, at
the end of the interval the exponential factor is
$\,\,-\frac{(E_{i}-M)^{2}}{2\sigma^{2}}=e^{-16\ln(10)}=10^{-16}$, which is
approximately equal to the machine precision for double accuracy used in our
calculations. We introduce a cut-off in the matrix size. We choose the
parameter $\sigma$ such that all eigenvalues in the region are good
approximations of the full eigenvalues. We will choose eigenvalues which
converge when the cut-off is increased. The precise definitions will be given
in the next section.

###  V.2 Eigenvalues convergence

In order to find the Gaussian states and their evolution, we look for some of
the eigenvalues and eigenvectors of the Hamiltonian operator. We considered a
cut-off: $n<1000,m<1000$ and calculated the matrix elements $H_{nm}$ using
equation (60). We calculated the integrals numerically using the DE rules for
infinite range integrals as they are described in [21]. This amounts to making
a change of variables

$q=e^{\pi\sinh(\tilde{q})},\quad p=\sinh(\pi\sinh(\tilde{p}))$ (112)

and performing the integrals in the range $]-4,4[$ using the open extended
trapezoidal rule (Trapz structure from [21]). We perform the integral over $q$
first and afterwards we perform the integral over $p$.

We look for convergence of the eigenvalues. Let us consider a family of
submatrices $H(k),k=1,\ldots,1000$, where each matrix $H(k)$ is obtained from
the Hamiltonian matrix $H$ by removing first $1000-k$ rows and first $1000-k$
columns. In particular $H(1000)=H$ and $k$ is the rank of the matrix $H(k)$.
Let us order the eigenvalues of $H(k)$ in the increasing order and let us
denote by $E_{i}(k),i\leq k$ the $i$-th eigenvalue of $H(k)$. We look for
convergence of the eigenvalues $E_{i}(k)$ as we increase $k$. The plot 3 shows
that eigenvalues $E_{i}$ for $700\leq i\leq 850$ converge. Each of the
eigenvalues stabilizes, which is reflected by the plateaux on the plot of
$E_{i}$ as a function of $k$.

Figure 3: The eigenvalues $E_{i}(k),i=700,710,\ldots,850$ of the submatrices
of $H(k)$. The eigenvalue $E_{i}(k)$ is the $i$-th eigenvalue of the sub-
matrix $H(k)$ (the eigenvalues are sorted in the increasing order). Since the
rank of $H(k)$ is equal $k$, it follows that $i\leq k$ and the plot of
$E_{i}(k)$ starts at $k=i$. For each $i$ we connected the values $E_{i}(k)$
with a line in order to guide the eye, i.e. each line represents one
eigenvalue $E_{i}$.

###  V.3 Evolution of the observables

We calculated the evolution of the expectation values of the operators
$\hat{q}$ and $\hat{p}$ in a coherent state peaked at energy $M=10^{-6}$ with
standard deviation $\sigma=10^{-7}$. This means that we considered a state:

$|\kern 1.29167pt\Psi(t)\kern
1.29167pt\rangle=\frac{1}{N}\sum_{i}e^{-\frac{(E_{i}-M)^{2}}{2\sigma^{2}}-\texttt{i}\,E_{i}\,t}\,|\kern
1.29167ptE_{i}\kern 1.29167pt\rangle$ (113)

and calculated the expectation values

$<\hat{q}>_{t}:=\langle\kern 1.29167pt\Psi(t)\kern 1.29167pt|\kern
1.29167pt\hat{q}\,\Psi(t)\kern 1.29167pt\rangle{\rm\ and\
}<\hat{p}>_{t}:=\langle\kern 1.29167pt\Psi(t)\kern 1.29167pt|\kern
1.29167pt\hat{p}\,\Psi(t)\kern 1.29167pt\rangle.$ (114)

The calculations are performed in the $e^{(1)}_{n}$ basis. After finding the
eigenvectors of $\hat{H}$, we calculated the components $\langle\kern
1.29167pte^{(1)}_{n}\kern 1.29167pt|\kern 1.29167pt\Psi(t)\kern
1.29167pt\rangle$. Afterwards, we used the formulas (109) and (100) for the
matrix elements of $\hat{q}$ and $\hat{p}$ in the $e^{(1)}_{n}$ basis.

We calculated $<\hat{q}>_{t}$ and $<\hat{p}>_{t}$ for $m=1.0,M=2\cdot
10^{-6},\sigma=10^{-7}$. The expectation value $<\hat{q}>_{t}$ is bounded from
below by approximately $0.006$. The bounce, connecting contracting and
expanding branches, occurs well above the horizon $q_{H}=2M^{2}=2\cdot
10^{-12}$ and at finite time as measured by an observer at spatial infinity.
We recall that in the corresponding classical case the shall falls onto the
horizon at an infinite time. The value $<\hat{q}>_{t}$ is also bounded from
above by approximately $0.01$, which is smaller than the classical maximal
size of the shell $q_{\rm max}=\frac{m^{4}}{8(m-M)^{2}}\approx\frac{1}{8}$. It
is clear that $<\hat{q}>_{t}$ must oscillate as illustrated by Fig. 4.

Since the classical variable $p$ is the measure of the angle between the
surfaces of constant time on both sides of the shell, it is reasonable to
expect that the expectation value $<\hat{p}>_{t}$ oscillates as a consequence
of the oscillation of $<\hat{q}>_{t}$. Indeed, such an effect is seen in Fig.
5.

Figure 4: Evolution of the expectation values of $\hat{q}$ in the coherent
state peaked at energy $M=10^{-6}$ with standard deviation $\sigma=10^{-7}$.
The expectation value is bounded from below by approximately $0.006$ ($0.037$
if the error bars are taken into account). This is well above the horizon
which in this case is at $q_{H}=2M^{2}=2\cdot 10^{-12}$. We interpret this as
a quantum bounce above the horizon. As a result, the evolution has oscillatory
character. Figure 5: Evolution of the expectation values of $\hat{p}$ in the
coherent state peaked at energy $M=10^{-6}$ with standard deviation
$\sigma=10^{-7}$. Let us notice that $<\hat{p}>_{t}=0$ when $<\hat{q}>_{t}$ is
extremal. Classically this happened only when $q$ is maximal.

###  V.4 Oscillatory character of the evolution

The oscillatory character of the evolution is related to our choice of the
Gaussian state. This can be shown by the following calculation:

$<\hat{q}>_{t}=\frac{1}{N^{2}}\sum_{i,j}e^{-\frac{(E_{i}-M)^{2}}{2\sigma^{2}}-\frac{(E_{j}-M)^{2}}{2\sigma^{2}}}e^{\texttt{i}(E_{i}-E_{j})t}\langle\kern
1.29167ptE_{i}\kern 1.29167pt|\hat{q}|\kern 1.29167ptE_{j}\kern
1.29167pt\rangle=\\\
=\frac{1}{N^{2}}\sum_{i,j}e^{-\frac{(E_{i}-M)^{2}}{2\sigma^{2}}-\frac{(E_{j}-M)^{2}}{2\sigma^{2}}}\cos((E_{i}-E_{j})t)\langle\kern
1.29167ptE_{i}\kern 1.29167pt|\hat{q}|\kern 1.29167ptE_{j}\kern
1.29167pt\rangle.$ (115)

In the formula above, we used the fact that the eigenvectors $|\kern
1.29167ptE_{i}\kern 1.29167pt\rangle$ have real components in the
$e^{(1)}_{n}$ basis and the fact $q_{nm}=q_{mn},\leavevmode\nobreak\
q_{nm}\in\mathbb{R}$. Let us recall that in our numerical calculations we
limited the sum to finite number of terms. Therefore the resulting function is
oscillatory as a finite sum of cosine functions.

Similar analysis applies to $<\hat{p}>_{t}$ giving trigonometric series.

We expect an oscillatory character if the spectrum is purely discrete. If the
spectrum turns out to be continuous, we expect that non-oscillatory functions
could be achieved. In that case the sum in (115) is replaced by an integral an
the resulting function may be non-oscillatory. A typical example of such
behaviour is a Fourier transform, which decomposes non-oscillatory function
into oscillatory modes. However, full spectral analysis is beyond the scope of
this paper.

##  VI Conclusions

We used the coherent state quantization technique to investigate a quantum
model of a massive shell. We developed the new mathematical tools: we found a
closed formula for the basis elements $e^{(1)}_{n}$, and for the matrix
elements of the operators $\hat{q}$ and $\hat{p}$ in this basis. We studied
the spectrum of the quantum Hamiltonian operator numerically in the case
$m=1.0$. We truncated the operator spectrum and observed that some of the
eigenvalues stabilize as we increase the truncation.

Using calculated eigenstates, we built a Gaussian state peaked at energy
$M=10^{-6}$. We investigated the evolution of the expectation values of
$\hat{q}$ and $\hat{p}$ in the Gaussian state and observed that the quantum
shell bounces well above the horizon. This leads to an oscillatory behaviour
of the system which is in contrast with the classical solution where the shell
collapses and reaches the Schwarzschild horizon $q_{H}=2M^{2}$ in infinite
time (as measured by an observer at infinity). The oscillatory character of
the evolution may be traced back to our construction of the Gaussian state.
More detailed analysis of this phenomenon needs a careful study of the
spectral properties of the Hamiltonian operator.

Our analysis of the quantum Hamiltonian is based on a series of numerical
experiments. The analytical analysis seems to be out of range. According to
these experiments it seems that this Hamiltonian, in the investigated region,
up to accuracy of our calculations (double accuracy) has a continuous
spectrum.

We realize that the research presented in this paper is not a complete
analysis. We developed the necessary tools and showed that this research
direction may lead to interesting results. In particular, the results
concerning oscillatory behaviour of the quantum system and the quantum origin
of a bounce above the horizon require independent cross-check. The integrals
involved in the expression of matrix elements of the Hamiltonian operator were
calculationally demanding. This resulted in rather limited size of the matrix
and small value of eigenvalues which can be reliably calculated. More
efficient numerical integration methods are needed in order to study the model
in more details. An advantage of our approach is that the quantum Hamiltonian
$\hat{H}$ we consider is positive definite, which is generally true and
independent on the chosen set of parameters.

The issue of the positivity of both classical and quantum Hamiltonians are of
special interest. At the classical level, it can be realized by restricting
the phase space. The quantum level is more demanding. In our integral
quantization method there is a one-to-one correspondence between the phase
space points and the space of coherent states. However, we cannot follow the
classical procedure of choosing the subspace of all coherent states to obtain
positivity of quantum Hamiltonian. The reason is that such procedure would
violate the mathematical consistency of the integral quantization and its
physical interpretation as deformation of POV measure which results in
problems with the probabilistic interpretation of this approach. Therefore, in
the integral defining the mapping of a classical observable into quantum
observable one cannot restrict the region of integration, but one can restrict
the integrant by taking $\leavevmode\nobreak\ \theta(H)H\leavevmode\nobreak\ $
instead of $H$. This way one obtains the quantum Hamiltonian $\hat{H}$ which
still acts on the full Hilbert space but it is positively defined because our
quantization procedure maps positive functions on the phase space into
positive operators. Similar logic has been used, for instance, in the paper
[22] within loop quantum gravity, where the dynamics is defined by a square
root of positive part of quantum gravitational scalar constraint operator.

To use this procedure consequently, the position and momentum quantum
observables have to be defined on the whole Hilbert space. They determine
physical meaning of the coherent states as quantum counterparts of the
configuration space points. Our procedure is consistent with treating the
operators $\hat{M}(Q)$, see Eq. (28), as the POV measure.

###### Acknowledgements.

We would like to thank Jacek Jezierski, Jerzy Kijowski, and Daniele Malafarina
for helpful discussions. This work was partially supported by the National
Science Centre, Poland grant No. 2018/28/C/ST9/00157.

##  Appendix A Orthonormal basis of the carrier space

The basis of the Hilbert space $L^{2}(\mathbb{R}_{+},d\nu(x))$, where
$d\nu(x):=dx/x$, is known to be [10]

$e^{(\alpha)}_{n}(x)=\sqrt{\frac{n!}{\Gamma(n+\alpha+1)}}\,e^{-x/2}x^{(1+\alpha)/2}\,L_{n}^{(\alpha)}(x),$
(116)

where $L_{n}^{(\alpha)}$ is the Laguerre function and $\alpha>-1$. One can
verify that
$\int_{0}^{\infty}e^{(\alpha)}_{n}(x)e^{(\alpha)}_{m}(x)d\nu(x)=\delta_{nm}$
so that $e^{(\alpha)}_{n}(x)$ is an orthonormal basis (for any fixed value of
the parameter $\alpha>-1$).

## References

  * [1] J. Kijowski, G. Magli, and D. Malafarina, “New derivation of the variational principle for the dynamikcs of a gravitational spherical shell”, Phys. Rev. D 74, 084017 (2006).
  * [2] C. Vaz, “Proper time quantization of a thin shell”, arXiv:2205.06867 [gr-qc].
  * [3] J. Jezierski and J. Kijowski, “Positivity of total energy in general relativity”, Phys. Rev. D 36, 1041 (1987).
  * [4] A. Góźdź, W. Piechocki and G. Plewa, “Quantum Belinski-Khalatnikov-Lifshitz scenario”, Eur. Phys. J. C 79, 45 (2019).
  * [5] A. Góźdź, W. Piechocki, and T. Schmitz, “Dependence of the affine coherent states quantization on the parametrization of the affine group”, Eur. Phys. J. Plus 136, 18 (2021).
  * [6] I. M. Gel′fand and M. A. Naïmark, “Unitary representations of the group of linear transformations of the straight line”, Dokl. Akad. Nauk. SSSR 55, 567 (1947).
  * [7] E. W. Aslaksen and J. R. Klauder, “Unitary Representations of the Affine Group”, J. Math. Phys. 9, 206 (1968).
  * [8] E. W. Aslaksen and J. R. Klauder, “Continuous Representation Theory Using Unitary Affine Group”, J. Math. Phys. 10, 2267 (1969).
  * [9] A. O. Barut and R. Ra̧czka, Theory of group representations and aplications (PWN, Warszawa, 1977).
  * [10] J. P. Gazeau and R. Murenzi, “Covariant affine integral quantization(s),” J. Math. Phys. 57, 052102 (2016).
  * [11] H. Bergeron and J. P. Gazeau, “Integral quantizations with two basic examples,” Annals Phys. 344, 43 (2014).
  * [12] M. Reed and B. Simon, Methods of Modern Mathematical Physics (San Diego, Academic Press, 1980), Vols I and II.
  * [13] G. B. Arfken, H. J. Weber, F. E. Harris, Mathematical Methods for Physicists: A Comprehensive Guide (Academic Press, Oxford, 2011).
  * [14] P. Busch, P.J. Lahti and P. Mittelstaedt, The Quantum Theory of Measurement, second Rev. Edition,Springer-Verlag, Berlin Heidelberg, 1996, ISBN 3-540-61355-2.
  * [15] A. Góźdź, A. Pȩdrak, and W. Piechocki, “Ascribing quantum system to Schwarzschild spacetime with naked singularity”, Class. Quantum Grav. 39, 145005 (2022).
  * [16] F. Faà di Bruno, “Note sur une nouvelle formule du calcul différentiel,” Quart. J. Math. , 1 (1855) pp. 359-360.
  * [17] Bell polynomial. Encyclopedia of Mathematics. URL: http://encyclopediaofmath.org/index.php?title=Bell_polynomial&oldid=46007
  * [18] E.T. Bell, “Exponential polynomials,” Ann. of Math. , 35 (1934) pp. 258-277.
  * [19] I. Lah, “Eine neue Art von Zahlen, ihre Eigenschaften und Anwendung in der mathematischen Statistik,” Mitteil. Math. Statist. , 7 (1955) pp. 203-216
  * [20] M. Abramowitz and I. A. Stegun, eds. Handbook of mathematical functions with formulas, graphs, and mathematical tables (US Government printing office, 1964), Vol. 55.
  * [21] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Brian, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 2007).
  * [22] T. Pawłowski and A. Ashtekar, “Positive cosmological constant in loop quantum cosmology”, Phys. Rev. D 85, 064001 (2012).


Observe that, if $\frac{1}{k}|Z-\bar{Z}|\leq\alpha_{k}$ is true, then we have

$\frac{1}{k}|f(Z)-f(\bar{Z})|\leq\frac{1}{k}L|Z-\bar{Z}|\leq L\alpha_{k}.$

Therefore we have

$\sup_{Z\in\mathcal{Z}_{k}}\mathbb{P}\left(\frac{1}{k}|f(Z)-f(\bar{Z})|\leq
L\alpha_{k}\right)\geq\sup_{Z\in\mathcal{Z}_{k}}\mathbb{P}\left(\frac{1}{k}|Z-\bar{Z}|\leq\alpha_{k}\right)\geq
1-q_{k},$

as claimed.

∎

We restate Lemma 32 below.

######  Lemma.

For each $k\geq 1$, let $\mathcal{I}_{k}$ be a set of tuples of nonnegative
random variables $(X,Y)\in\mathcal{I}_{k}$, with deterministic fluid
approximations $(\bar{X},\bar{Y})$, which satisfy the asymptotic concentration
property (98). Let $(M,R,p)$ be any constants satisfying $M\leq\gamma k$,
$R\leq\gamma k$, and $p\in[0,1]$. Let $G$ be the number of dispatches and $U$
the number of remaining drivers, when the single-destination dispatch
subroutine (1) is used to allocate $R-Y$ dispatch requests to $M-X$ drivers,
assuming that, except for a subset of size at most $k\delta_{k}$ drivers,
drivers use an acceptance probability within $\epsilon_{k}$ of $p$. Let
$\bar{G}$ and $\bar{U}$ be the output of the fluid subroutine when $R-\bar{Y}$
riders are allocated to $M-\bar{X}$ drivers. Then $G$ converges asymptotically
to $\bar{G}$ and $U$ converges asymptotically to $\bar{U}$, in the sense that
the following equation holds for all $k\geq 1$:

$\displaystyle\sup_{M\leq\gamma k,\ R\leq\gamma k,\ p\in[0,1],\
(X,Y)\in\mathcal{I}_{k}}\mathbb{P}\left(\frac{1}{k}|G-\bar{G}|\leq\beta_{k}\right)\geq
1-p_{k}$ (108) $\displaystyle\sup_{M\leq\gamma k,\ R\leq\gamma k,\ p\in[0,1],\
(X,Y)\in\mathcal{I}_{k}}\mathbb{P}\left(\frac{1}{k}|U-\bar{U}|\leq\beta_{k}\right)\geq
1-p_{k},$ (109)

where $(\beta_{k}:k\geq 1)$ and $(p_{k}:k\geq 1)$ are nonnegative sequences
which converge to $0$ as $k\to\infty$.

###### Proof.

Consider the single-destination dispatch subroutine when the population-size
parameter is $k\geq 1$, with $M\leq\gamma k$ drivers, $R\leq\gamma k$ riders,
$p\in[0,1]$ common acceptance probability, and let $(X,Y)\in\mathcal{I}_{k}$.
Let $G^{\prime}$ and $U^{\prime}$ be the output of the fluid matching process
with $R-Y$ riders and $M-X$ drivers. By Lemma 31 we have that $G$ and $U$
converge asymptotically to $G^{\prime}$ and $U^{\prime}$, in the sense that
there exists $(\alpha^{\prime}_{k}:k\geq 1)$, $(q^{\prime}_{k}:k\geq 1)$, both
converging to $0$ as $k\to\infty$, such that the following holds for every
$k$:

$\displaystyle\sup_{M\leq\gamma k,\ R\leq\gamma k,\ p\in[0,1],\
(X,Y)\in\mathcal{I}_{k}}\mathbb{P}\left(\frac{1}{k}|G-G^{\prime}|\leq\alpha^{\prime}_{k}\right)\geq
1-q^{\prime}_{k}$ $\displaystyle\sup_{M\leq\gamma k,\ R\leq\gamma k,\
p\in[0,1],\
(X,Y)\in\mathcal{I}_{k}}\mathbb{P}\left(\frac{1}{k}|U-U^{\prime}|\leq\alpha^{\prime}_{k}\right)\geq
1-q^{\prime}_{k}.$

We now claim that $G^{\prime}$ and $U^{\prime}$ converge asymptotically to
$\bar{G}$ and $\bar{U}$. Recall the output of the fluid single-destination
dispatch subroutine is defined (Definition 4) in terms of deterministic
functions $\bar{G}(\cdot)$ and $\bar{U}(\cdot)$, so we have

$\displaystyle G^{\prime}=\bar{G}(R-Y,M-X,p)$
$\displaystyle\bar{G}=\bar{G}(R-\bar{Y},M-\bar{X},p)$ $\displaystyle
U^{\prime}=\bar{U}(R-Y,M-X,p)$
$\displaystyle\bar{U}=\bar{U}(R-\bar{Y},M-\bar{X},p).$

Observe that the functions $\bar{G}(\cdot)$ and $\bar{U}(\cdot)$ are both
Lipschitz continuous, so by Lemma 34 we have $G^{\prime}$ converges
asymptotically to $\bar{G}$, i.e. the following equation holds for every
$k\geq 1$

$\sup_{M\leq\gamma k,\ R\leq\gamma k,\ p\in[0,1],\
(X,Y)\in\mathcal{I}_{k}}\mathbb{P}\left(\frac{1}{k}|G^{\prime}-\bar{G}|\leq
L\alpha_{k}\right)\geq 1-q_{k},$

where $\alpha_{k}$ and $q_{k}$ are the error term and probability term from
the convergence of $(X,Y)$ to $(\bar{X},\bar{Y})$ (see equation (98)), and $L$
is the Lipschitz constant for the function $\bar{G}(\cdot)$. Taking
$\beta_{k}=\alpha^{\prime}_{k}+L\alpha_{k}$ and $p_{k}=q^{\prime}_{k}+q_{k}$
finishes the proof. ∎

###  E.5 Proof of Lemma 33

We restate Lemma 33 below.

######  Lemma.

Suppose that the total number of relocation trips converges to the fluid
volume of relocation trips, as $k\to\infty$, for all admissible driver states.
That is, assume there exists sequences $(\alpha_{k}:k\geq 1)$ and
$(q_{k}:k\geq 1)$ such that
$\lim_{k\to\infty}\alpha_{k}=\lim_{k\to\infty}q_{k}=0$, for which the
following inequality holds for every $k$ and every $\ell$

$\sup_{\mathcal{M}\in\mathcal{D}_{k}}\mathbb{P}\left(\frac{1}{k}|H_{\ell}-\bar{H}_{\ell}|\leq\alpha_{k}\right)\geq
1-q_{k}.$

Then the relocation trip volumes along each individual route converge to their
fluid approximations, i.e. there exist sequences $(\beta_{k}:k\geq 1)$ and
$(p_{k}:k\geq 1)$ such that
$\lim_{k\to\infty}\beta_{k}=\lim_{k\to\infty}p_{k}=0$, for which the following
inequality holds for every $k$:

$\sup_{\mathcal{M}\in\mathcal{D}_{k}}\mathbb{P}\left(\frac{1}{k}\|\bm{H}-\bar{\bm{H}}\|_{1}\leq\beta_{k}\right)\geq
1-p_{k}.$

###### Proof.

It suffices to prove the following is true for each route $(\ell,d)$:

$\sup_{\mathcal{M}\in\mathcal{D}_{k}}\mathbb{P}\left(\frac{1}{k}\|H_{(\ell,d)}-\bar{H}_{(\ell,d)}\|_{1}\leq\beta_{k}\right)\geq
1-p_{k}.$

If the above inequality holds for each route $(\ell,d)$, then the claimed
inequality follows from a union bound over all locations.

Fix $k\geq 1$, an admissible driver state $\mathcal{M}\in\mathcal{D}_{k}$, and
a location $\ell$. Let $\mathcal{M}_{\ell}^{R}\subseteq\mathcal{M}$ be the
stochastic subset of drivers who serve a relocation trip. By definition, we
have $|\mathcal{M}_{\ell}^{R}|$ is the random variable $H_{\ell}$, and by
assumption we have that $H_{\ell}$ converges to a deterministic fluid
approximation $\bar{H}_{\ell}$.

The probability any individual driver falls in this subset depends on the
exact threshold vector that the driver has selected, as well as their realized
add-passenger disutility. In particular, if two drivers use the exact same
threshold vector, they have the same probability of going non-dispatched. If
two drivers use approximately the same threshold, and we condition on the
event that their sampled add-passenger disutilities are bounded away from the
region where the different thresholds would lead to different decisions, then
drivers still have the same probability of going non-dispatched.

Let $\bm{x}_{\ell}$ be the common disutility threshold vector which is
approximately used by approximately all drivers at $\ell$. Define

$\mathcal{X}(\epsilon_{k},\bm{x}_{\ell})=[0,C]\setminus\left(\bigcup_{d\in\mathcal{L}}[x_{(\ell,d)}-\epsilon_{k},x_{(\ell,d)}+\epsilon_{k}]\right)$

to be the subset of feasible disutility thresholds $[0,C]$ where a small band
$[x_{(\ell,d)}-\epsilon_{k},x_{(\ell,d)}+\epsilon_{k}]$ centered at each
threshold $x_{(\ell,d)}$ is removed.

Define

$\mathcal{M}_{\ell}^{\prime}=\mathcal{M}_{\ell}(\epsilon_{k},\bm{x}_{\ell})\cap\\{i\in\mathcal{M}_{\ell}:X_{i}\in\mathcal{X}(\epsilon_{k},\bm{x}_{\ell})\\}$

to be the subset of drivers who approximately use the threshold
$\bm{x}_{\ell}$ and whose sampled disutilities lie in
$\mathcal{X}(\epsilon_{k},\bm{x}_{\ell})$. Notice that the cardinality
$|\mathcal{M}_{\ell}^{\prime}|$ has Binomial distribution with parameters
$|\mathcal{M}_{\ell}(\epsilon_{k},\bm{x}_{\ell})|$ and
$1-2|\mathcal{L}|\epsilon_{k}$, so $|\mathcal{M}_{\ell}^{\prime}|$ converges
to $|\mathcal{M}_{\ell}|$ as $k\to\infty$.

Define

$H^{\prime}_{(\ell,d)}=\left|\\{i\in\mathcal{M}_{\ell}:r_{i}=d\\}\cap\mathcal{M}_{\ell}^{\prime}\cap\mathcal{M}_{\ell}^{R}\right|.$

Notice that the number of trips from $\ell$ to $d$ is bounded by

$|\mathcal{M}_{\ell}\setminus\mathcal{M}_{\ell}^{\prime}|+H^{\prime}_{(\ell,d)}\geq
H_{(\ell,d)}\geq H^{\prime}_{(\ell,d)}.$

Finally, notice that the distribution of $H^{\prime}_{(\ell,d)}$ is equivalent
to sampling $|\mathcal{M}_{\ell}^{R}|$ balls, without replacement, from a bag
with $|\mathcal{M}^{\prime}_{\ell}|$ balls, where each ball is associated with
a destination, and counting how many balls are associated with the destination
$d$. Concentration inequalities for sampling without replacement show that
$H^{\prime}_{(\ell,d)}$ converges asymptotically to $\bar{H}_{(\ell,d)},$ and
this is sufficient to prove that $H_{(\ell,d)}$ converges to
$\bar{H}_{(\ell,d)}$ asymptotically, because the difference between
$H^{\prime}_{(\ell,d)}$ and $H_{(\ell,d)}$ vanishes asymptotically.

∎

##  Appendix F General Convex Analysis Properties

In this section we obtain useful convex analysis properties. We change
notation from the rest of the paper, and consider the following generic convex
optimization problem

$\inf\left(f(x)\mid g_{i}(x)\leq 0,i=1,2,\dots,m,x\in\mathbb{R}^{n}\right),$
(110)

where $f,g_{1},\dots,g_{m}$ are convex functions from $\mathbb{R}^{n}$ to
$\mathbb{R}$. We assume our convex program (110) satisfies the conditions
described in Assumption 5.

######  Assumption 5.

Assume the following conditions hold:

  1. 1.

$f,g_{1},\dots,g_{m}$ are continuously differentiable at every point
$x\in\mathbb{R}^{n}$.

  2. 2.

The feasible region $\\{x\mid g_{i}(x)\leq
0,i=1,2,\dots,m,x\in\mathbb{R}^{n}\\}$ is a bounded compact set. In
particular, there is a constant $\gamma$ such that the feasible region lies in
$\bar{B}(0,\gamma)$, i.e. the closed ball centered at $0$ with radius
$\gamma$.

  3. 3.

The gradients $\nabla f(x),\nabla g_{1}(x),\dots,\nabla g_{m}(x)$ have norm
smaller than some constant $C$ for all $x\in\mathbb{R}^{n}$.

The Lagrangian associated with the optimization problem (110) is the function
$L:\mathbb{R}^{n}\times\mathbb{R}^{m}_{+}\to\mathbb{R}\cup\\{+\infty\\}$
defined by

$L(x;\lambda)=f(x)+\lambda^{T}g(x),$ (111)

where $g(x)$ is the vector in $\mathbb{R}^{m}$ with $g_{i}(x)$ as its $i$th
component. Note the min-max theorem states the following relads :

$\inf_{x\in\mathbb{R}^{n}}\sup_{\lambda\in\mathbb{R}^{m}_{+}}L(x;\lambda)=\sup_{\lambda\in\mathbb{R}^{m}_{+}}\inf_{x\in\mathbb{R}^{n}}L(x;\lambda).$
(112)

######  Definition 8.

Let $x\in\mathbb{R}^{n}$. A nonnegative vector $\lambda\in\mathbb{R}^{m}_{+}$
is said to be a Lagrange multiplier vector for $x$ if it satisfies the
following conditions:

  * •

Complementary slackness: $\lambda_{i}g_{i}(x)=0$ holds for every
$i=1,2,\dots,m$.

  * •

Stationarity: $\nabla_{x}L(x;\lambda)=0$. In other words, $x$ is a global
minimizer of $L(\cdot;\lambda)$.

We say $\lambda$ is an $\epsilon$-approximate Lagrange multiplier vector for
$x$ if it satisfies the following conditions:

  * •

Approximate complementary slackness: $|\lambda_{i}g_{i}(x)|\leq\epsilon$ holds
for every $i=1,2,\dots,m$.

  * •

Approximate stationarity: $\|\nabla_{x}L(x;\lambda)\|_{2}\leq\epsilon$.

For the next Lemma, we consider an optimal solution $x^{*}$ and a Lagrange
multiplier vector $\lambda^{*}$ for $x^{*}$. Part of our proof is concerned
with the function mapping $x\in\mathbb{R}^{n}$ to the gradient of the
Lagrangian $\nabla_{x}L(x;\lambda^{*})$. In particular, we care about how
large the gradient can vary when evaluated at two points that are close to one
another. To reason about this maximum perturbation effect, for $\delta>0$
define

$\kappa(\delta)=\max\left(\|\nabla_{x}L(x;\lambda^{*})-\nabla_{x}L(y;\lambda^{*})\|_{2}\mid
x,y\in\bar{B}(0,\gamma),\|x-y\|_{2}\leq\delta\right)$ (113)

which gives the maximum norm of the difference between the gradients of any
two points in $\bar{B}(0,\gamma)$ whose distance from each other is at most
$\delta$. From Assumption 5 we know the gradient $\nabla_{x}L(x;\lambda^{*})$
is continuous in $x$, and we know that every continuous function over a
compact set is uniformly continuous, so it follows that the gradient function,
restricted to the closed ball $\bar{B}(0,\gamma)$, is uniformly continuous.
Therefore, the maximum perturbation $\kappa(\delta)$ goes to $0$ as $\delta\to
0$.

######  Lemma 35.

Let $\bar{x}$ be a feasible solution for (110) and let
$\bar{\lambda}\in\mathbb{R}^{m}_{+}$. If $\bar{\lambda}$ is a Lagrange
multiplier vector for $\bar{x}$ then $\bar{x}$ is an optimal solution for
(110). If $\bar{\lambda}$ is an $\epsilon$-approximate Lagrange multiplier
vector for $\bar{x}$ then $\bar{x}$ is an $\epsilon^{\prime}$-optimal solution
for (110), where $\epsilon^{\prime}=\epsilon(n+2\gamma)$.

###### Proof.

We focus on the case where $\bar{\lambda}$ is an $\epsilon$-approximate
Lagrange multiplier vector for $\bar{x}$. Let $x^{*}$ and $\lambda^{*}$ be
optimal primal and dual variables. Observe the following chain of
inequalities:

$f(x^{*})=L(x^{*};\lambda^{*})\geq L(x^{*};\bar{\lambda})\geq
L(\bar{x};\bar{\lambda})+\nabla_{x}L(\bar{x};\bar{\lambda})^{T}(x^{*}-\bar{x}),$

where the first bound follows from optimality of $\lambda^{*}$ and the second
line follows by convexity of $L(\cdot;\bar{\lambda})$. Therefore we have the
upper bound

$f(\bar{x})\leq
f(x^{*})+\|\nabla_{x}L(\bar{x};\bar{\lambda})\|_{2}\|(x^{*}-\bar{x})\|_{2}+|\bar{\lambda}^{T}g(\bar{x})|.$

By approximate complementary slackness we have
$|\bar{\lambda}^{T}g(\bar{x})|\leq n\epsilon$, by approximate stationarity we
have $\|\nabla_{x}L(\bar{x};\bar{\lambda})\|_{2}\leq\epsilon$, and by
feasibility of $x^{*}$ and $\bar{x}$ we have $\|(x^{*}-\bar{x})\|_{2}\leq
2\gamma$. Therefore we obtain the bound

$f(\bar{x})\leq f(x^{*})+\epsilon(n+2\gamma)=f(x^{*})+\epsilon^{\prime}$

as claimed. ∎

######  Lemma 36.

Let $x^{*}$ be an optimal solution for (110) and let $\lambda^{*}$ be a
Lagrange multiplier vector for $x^{*}$. Let $\bar{x}$ be a feasible solution
for (110). If $\bar{x}$ is an optimal solution then $\lambda^{*}$ is a
Lagrange multiplier vector for $\bar{x}$. If $\bar{x}$ is an
$\epsilon$-optimal solution, then $\lambda^{*}$ is an
$\epsilon^{\prime}$-approximate Lagrange multiplier vector, where
$\epsilon^{\prime}=\sqrt{\sqrt{\epsilon}+C\kappa(C\sqrt{\epsilon})}$.

###### Proof.

Observe the following chain of inequalities:

$\displaystyle f(x^{*})$ $\displaystyle=$ $\displaystyle L(x^{*};\lambda^{*})$
$\displaystyle=$ $\displaystyle\min_{x\in\mathbb{R}^{n}}L(x;\lambda^{*})$
$\displaystyle\leq$ $\displaystyle L(\bar{x};\lambda^{*})$ $\displaystyle=$
$\displaystyle f(\bar{x})+\lambda^{*T}g(\bar{x})$ $\displaystyle\leq$
$\displaystyle f(\bar{x})$ (115)

The first two lines follow from the complementary slackness and stationarity
conditions which hold between $x^{*}$ and $\lambda^{*}$. Line (115) follows
from feasibility of $\bar{x}$ and nonnegativity of $\lambda^{*}$.

We first consider the case where $\bar{x}$ is an optimal solution. In this
case, we have $f(\bar{x})=f(x^{*})$ so every line in the above chain of
inequalities holds with equality. In particular, line (F) holding with
equality shows that $\bar{x}$ is a global minimizer of $L(\cdot;\lambda^{*})$,
so the stationarity condition $\nabla_{x}L(\bar{x};\lambda^{*})$ holds, and
line (115) holding with equality shows that complementary slackness holds
between $\bar{x}$ and $\lambda^{*}$. Therefore $\lambda^{*}$ is a Lagrange
multiplier vector for $\bar{x}$.

Next, we consider the case wehre $\bar{x}$ is an $\epsilon$-optimal solution.
In this case, we have $f(\bar{x})\leq f(x^{*})+\epsilon$, so $\epsilon$ is an
upper bound on the difference between any two consecutive terms in the above
chain of inequalities. Applying this upper bound to line (115) we conclude
that $|\lambda^{*T}g(\bar{x})|\leq\epsilon$. Since every term
$\lambda^{*}_{i}g_{i}(\bar{x})$ is nonpositive, it follows that
$|\lambda^{*}_{i}g_{i}(\bar{x})|\leq|\lambda^{*T}g(\bar{x})|\leq\epsilon\leq\epsilon^{\prime}$,
so approximate complementary slackness holds between $\lambda^{*}$ and
$\bar{x}$.

We show in Lemma 37 that $\bar{x}$ satisfies approximate stationarity with
respect to $\lambda^{*}$. We use the result in Lemma 37 by taking $h(x)$ to
mean $L(x;\lambda^{*})$. The result of Lemma 37 then provides the bound

$\|\nabla_{x}L(\bar{x};\lambda^{*})\|_{2}^{2}\leq\sqrt{\epsilon}+C\kappa(C\sqrt{\epsilon}).$

Therefore, $\|\nabla_{x}L(\bar{x};\lambda^{*})\|_{2}\leq\epsilon^{\prime}$,
establishing that $\lambda^{*}$ is an $\epsilon^{\prime}$-approximate Lagrange
multiplier vector for $\bar{x}$. ∎

######  Lemma 37.

Let $h:\mathbb{R}^{n}\to\mathbb{R}$ be a continuously differentiable convex
function. Let $x^{*}$ be a global minimizer of $h$, let $\bar{x}$ satisfy
$h(\bar{x})<h(x^{*})+\epsilon$, and assume that $\bar{x}$ and $x^{*}$ both lie
in $\bar{B}(0,\gamma)$, i.e. the closed ball centered at $0$ with radius
$\gamma$. For $\delta>0$, define

$\kappa(\delta)=\max\left(\|\nabla h(x)-\nabla h(y)\|_{2}\mid
x,y\in\bar{B}(0,\gamma),\|x-y\|_{2}\leq\delta\right)$

to be the maximum norm of the difference between the gradients of any two
points in $\bar{B}(0,\gamma)$ whose distance from each other is at most
$\delta$. Assume $C$ is a constant upper bound on $\|\nabla h(x)\|_{2}$. Then

$\|\nabla h(\bar{x})\|_{2}^{2}\leq\sqrt{\epsilon}+C\kappa(C\sqrt{\epsilon}).$
(116)

###### Proof.

Define

$x^{\prime}=\bar{x}-t\nabla h(\bar{x}).$

for some $t>0$. From convexity of $h$ we have

$h(x^{\prime})+\langle\nabla h(x^{\prime}),\bar{x}-x^{\prime}\rangle\leq
h(\bar{x}).$

Rearranginge the above, and using the definition of $x^{\prime}$,

$t\langle\nabla h(x^{\prime}),\nabla h(\bar{x})\rangle\leq
h(\bar{x})-h(x^{\prime})\leq\epsilon.$

Now consider the following chain of inequalities:

$\displaystyle\|\nabla h(\bar{x})\|_{2}^{2}$ $\displaystyle=$
$\displaystyle\langle\nabla h(\bar{x}),\nabla h(\bar{x})\rangle$
$\displaystyle=$ $\displaystyle\langle\nabla h(x^{\prime}),\nabla
h(\bar{x})\rangle+\langle\nabla h(\bar{x})-\nabla h(x^{\prime}),\nabla
h(\bar{x})\rangle$ $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{t}+\langle\nabla h(\bar{x})-\nabla
h(x^{\prime}),\nabla h(\bar{x})\rangle$ $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{t}+\|\nabla h(\bar{x})-\nabla
h(x^{\prime})\|_{2}\|\nabla h(\bar{x})\|_{2}$ $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{t}+\|\nabla h(\bar{x})-\nabla
h(x^{\prime})\|_{2}C$ $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{t}+C\kappa(\|\bar{x}-x^{\prime}\|_{2})$
$\displaystyle\leq$ $\displaystyle\frac{\epsilon}{t}+C\kappa(\|t\nabla
h(\bar{x})\|_{2})$ $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{t}+C\kappa(Ct).$

Setting $t=\sqrt{\epsilon}$, we obtain the claimed bound (116). ∎