Classical and quantum billiard inside the square with gravitational field

During the last decades dynamical systems, especially billiard systems, have widely been studied in various context such as the theory of classical periodic orbits [8, 10] in polygonal systems, quantum versions and visualizations [6] or applications to number theory [3]. Due to simplicity those billiard systems are often studied in a force free scenario where the particles motion is determined by its initial conditions as well as the law of reflection at the boundary. The next interesting systems are those involving a constant force along one direction which e.g. can be induced by a gravitational field. Recent research on this topic has e.g. classically been performed for the case of a parabolic boundary in [4] or in the quantum version either analytically in the simplest case [2] or numerically in [1].

In this research we investigate the classical and quantum behavior of a particle inside a square box which is known to be integrable.

The structure of the paper is as follows: In section 2 we consider the classical motion of a particle inside a square with gravitational field along the $y-$direction and derive the conditions for periodic orbits. In the following section we study the classical probability density for this setup and derive the expected position as well as the standard deviation of the particle including a study of limiting cases for small and large vertical energies. Section 4 will be focusing on the quantum mechanical description resulting in a derivation of the associated wave functions and quantized energy eigenvalues. The paper concludes in the quantum mechanical study of the probability density and the associated expectation values for the position, including a comparison to the classical results.

It turns out that periodic orbits are completely determined by the trajectory behavior along the boundaries. For a particle moving inside a square of side length $L$ under the influence of gravity along the $y-$axis the potential is given by

$$V(x,y)=\begin{cases}mgy&\mbox{for}~{}(x,y)\in[0;L]^{2},\\ \infty&\mbox{otherwise}.\end{cases}$$

(1)

The motion of the particle is totally characterized by the initial starting position $x_{0}$, the absolute value of the velocity $v$ (or equivalently its total energy $E$) at $y=0$ and the angle $\varphi$ under which the particle ejects into the domain (see Figure 1).

The equations of motion for time $t$ before the first scattering along the boundary are given by

	$\displaystyle x(t)$	$\displaystyle=v\cdot\cos(\varphi)\cdot t+x_{0},$		(2)
	$\displaystyle y(t)$	$\displaystyle=-\frac{1}{2}g\cdot t^{2}+v\cdot\sin(\varphi)\cdot t.$		(3)

It turns out that rewriting $y$ in terms of $x$ and replacing the velocity dependence by the associated total energy of the system $E=\frac{1}{2}mv^{2}$ is an advantage. As a result, one obtains the particle height in terms of its distance along the $x-$direction via

$$y(x)=-\frac{mg}{4E\cdot\cos^{2}(\varphi)}\cdot(x-x_{0})^{2}+\tan(\varphi)\cdot(x-x_{0}).$$

(4)

Trajectories thus correspond to parabolas. Hitting the left or right boundary simply corresponds to inverting the $x-$direction of the trajectory. The particle motion can thus be alternatively described by unfolding the square [3, 8, 10] and consider the alternating parabola trajectories in the unfolded picture (see Figure 2).

For energies $E<\frac{mgL}{\sin^{2}(\varphi)}$, which simply means that the particle energy in the $y-$direction is not enough to encounter the upper boundary and thus only propagating in domains $0\leq y<L$, periodic orbits thus arise if the difference between two consecutive zeros of $y(x)$, i.e. $\Delta x(y=0)$, satisfies

$$p\cdot\Delta x(0)=q\cdot 2L~{}~{}~{}~{}~{}p,q\in\mathbb{N}^{+}.$$

(5)

The factor 2 represents the fact that the domain has to be unfolded twice to return to its initial orientation (see Figure 2).

If the particles energy along the $y-$direction is large enough to obtain (in theory) heights ${y\geq L}$, the particle bounces at the upper boundary. Since in this case the $y-$component of the velocity is inverted ($v_{y}\rightarrow-v_{y}$) the effective distance covered by two consecutive bounces at $y=0$ with one intermediate bounce at $y=L$ (see Figure 3) is given by $\Delta x(0)-\Delta x(L)$, where

$$\Delta x(y)=\frac{2E}{mg}\cdot\sin(2\varphi)\cdot\sqrt{1-\frac{mgy}{E\cdot\sin^{2}(\varphi)}}.$$

(6)

Applying the unfolding procedure we can see that the particle describes periodic orbits iff the difference $\Delta x(0)-\Delta x(L)$ satisfies

$$[\Delta x(0)-\Delta x(L)]\cdot p=q\cdot 2L~{}~{}~{}~{}p,q\in\mathbb{N},$$

(7)

i.e. periodic orbits in the square geometry with gravity are completely determined by the behavior of the particle along the boundaries at $y=0$ and $y=L$. Note that (7) also covers the first case displayed in Eq. (5) with energies along the $y-$direction below $mgL$. In that scenario, it is implicitly understood that $\Delta x(L)=0$ if $E\cdot\sin^{2}(\varphi)\leq mgL$. For the particle trajectories not hitting a corner of the square, in which case the law of reflection is not defined, there are two further conditions - namely

$$x_{0}+[\Delta x(0)-\Delta x(L)]\cdot(n+\delta)\neq m\cdot L,$$

(8)

with $n,m\in\mathbb{N}$ and $\delta=0$ for hitting the corners $A$ or $B$ and $\delta=\frac{1}{2}$ for corners $C$ or $D$.

In this section we are first considering the probability density of a particle moving inside the square including gravity along the $y-$direction. The maximal theoretical height (neglecting the boundary) a particle might reach is given in terms of the total energy $E$, its mass $m$ and the initial angle $\varphi$ by

$$h_{max}=\frac{E\cdot\sin^{2}(\varphi)}{mg}.$$

(9)

Note that in the following simply abbreviate $h_{max}=h$. From here we start off by considering the classical probability density $\varrho(x,y)$ which has been discussed [5, 7] in various other physical systems e.g. free particles or the harmonic oscillator. Due to the fact that the motion in the $x-$direction is homogeneous, there is no dependence on $x$ in the probability density, i.e. $\varrho(x,y)=\varrho(y)$. Classically, the probability of finding the particle in the vicinity $dx\cdot dy$ is proportional to the differential time $dt$ of the particle in this vicinity, i.e.

$$\varrho(y)dxdy\sim dt.$$

(10)

The last equation is equivalent to

$$\varrho(y)dxdy\sim\frac{dy}{\frac{dy}{dt}}=\frac{dy}{v_{y}(y)}.$$

(11)

By conservation of energy, it is easy to verify that the velocity $v_{y}$ in $y-$direction depending on the energy $E$ and potential energy $V(y)=mgy$ of the particle at a height $y$ is proportional to

$$v_{y}(y)\sim\sqrt{E\cdot\sin^{2}(\varphi)-mgy}.$$

(12)

A direct calculation shows that the classical probability density of the particle written in terms of the maximal theoretical height $h$ is given by

$$\varrho(y)=\frac{N}{\sqrt{h-y}}\cdot\Theta\left(h-y\right),$$

(13)

where the Heaviside Theta-function is defined by

$$\Theta(x)=\begin{cases}0&\mbox{for}~{}x<0,\\ 1&\mbox{for}~{}x\geq 0,\end{cases}$$

(14)

and insures that the $y-$integration runs from zero to the maximal allowed height $h$ but not further than the boundary of the domain $L$. The correct normalization constant $N$ can, as usual, be obtained by the condition

$$\int_{0}^{L}\int_{0}^{L}\varrho(y)dydx=1,$$

(15)

and takes the form

$$N=\begin{cases}\frac{1}{2L\cdot\sqrt{h}}&\mbox{for}~{}h<L,\\ \frac{1}{2L\cdot\left(\sqrt{h}-\sqrt{h-L}\right)}&\mbox{for}~{}h\geq L.\end{cases}$$

(16)

Figure 4 displays some cases for the obtained probability density along the $y-$direction for
$L=1$ and various vertical energies. Obviously, the probability density approaches $1/L^{2}$ for large energies (green curve) due to the reflection along the upper boundary and displays an asymptotic behavior (orange curve) in the low energy limit since in this regime $y<L$.

Applying Eq. (13) with associated normalization Eq. (16) it is now straight forward to calculate the classical expected particle positions $\langle\vec{r}\rangle_{cl}$ inside the box as a function of $h$ (i.e. corresponding energy) via

$$\langle\vec{r}\rangle_{cl}=\int_{0}^{L}\int_{0}^{L}\vec{r}\cdot\varrho(y)dxdy.$$

(17)

Note that the same expression can later on be applied for the quantum case study where the probability density is given by the absolute square of the wave function. Similarly, one can calculate the expectation values $\langle x^{2}\rangle_{cl}$ and $\langle y^{2}\rangle_{cl}$ which are needed for the calculation of the particles uncertainty (or standard deviation ${\Delta x_{cl}=\sqrt{\langle x^{2}\rangle_{cl}-\langle x\rangle_{cl}^{2}}}$ and analogous for $\Delta y_{cl}$.

We will begin with a short discussion of the $x-$direction. Since in this direction there is no force acting on the particle, the probability density and associated expectation value calculations can easily be carried out, yielding

$$\langle x\rangle_{cl}=\frac{L}{2},$$

(18)

and

$$\langle x^{2}\rangle_{cl}=\frac{L^{2}}{3}.$$

(19)

As a direct result, the classical uncertainty of the particle along the $x-$direction is given by

$$\Delta x_{cl}=\frac{L}{2\sqrt{3}}.$$

(20)

The results for the $y-$direction can also be performed, resulting in the expressions (again in terms the Heaviside function)

$$\langle y\rangle_{cl}=\frac{2}{3}h-\frac{L}{3\cdot\left(\sqrt{\frac{h}{h-L}}-1\right)}\cdot\Theta(h-L),$$

(21)

and

$$\langle y^{2}\rangle_{cl}=\frac{8h^{2}}{15}-\frac{4hL+3L^{2}}{15\cdot\left(\sqrt{\frac{h}{h-L}}-1\right)}\cdot\Theta(h-L).$$

(22)

Remarkably, the result for $\langle y\rangle_{cl}$ demonstrates that the maximal expected height of the particle is given by $\frac{2}{3}h$ and is accomplished if ${E\cdot\sin^{2}(\varphi)=mgL}$, i.e. if the vertical energy of the particle matches the potential energy at the upper boundary. For energies above $mgL$ the expectation value $\langle y\rangle_{cl}$ approaches $L/2$ which is the same as $\langle x\rangle_{cl}$ and is clear from the fact that in this limit the gravitational force is merely neglectable resulting basically in the movement of a free particle inside a square box. For the classical standard deviation this results in

$$\Delta y_{cl}=\frac{2\sqrt{15}}{15}\cdot h\cdot\sqrt{1+J(L,h)\cdot\Theta(h-L)},$$

(23)

where, in order to write the formula in a compact form, we defined the new function $J(L,h)$ corresponding to the correction terms for the standard deviation induced by the upper boundary compared to the system without the boundary by

$$J=\frac{L\left(8h\left(\sqrt{\frac{h}{h-L}}-1\right)+L\left(4-9\sqrt{\frac{h}{h-L}}\right)\right)}{4h^{2}\left(\sqrt{\frac{h}{h-L}}-1\right)^{2}}.$$

(24)

Again, as for the expected particle position in $y-$direction in the large vertical energy limit, the standard deviation approaches the same value as in $x-$direction, i.e.

$$\lim_{h\rightarrow\infty}\Delta y_{cl}=\Delta x_{cl}=\frac{L}{2\sqrt{3}}.$$

(25)

A graphical representation for $\langle y\rangle_{cl}$ (blue curve) and the associated uncertainty $\Delta y_{cl}$ (orange domain) as a function of $h$ and the special case of a unit square with $L=1$ is displayed in Figure 5.

The Schrödinger equation for a quantum particle of mass $m$ inside a square of length $L$ and potential given by Eq. (1) reads

$$\left[-\frac{\hbar^{2}}{2m}\cdot\Delta_{2}+mgy\right]\cdot\Psi(x,y)=E\cdot\Psi(x,y),$$

(26)

where $\Delta_{2}=\partial_{x}^{2}+\partial_{y}^{2}$ is the two dimensional Laplace-operator. Referring to the classical system, where the $x-$ and $y-$motion are independent from each other the wave function can be represented by the product of two separate functions for each variable, i.e. $\Psi(x,y)=X(x)\cdot Y(y)$. With this ansatz the Schrödinger equation basically falls apart into two independent equations for $X(x)$ and $Y(y)$ separately, namely

	$\displaystyle-\frac{\hbar^{2}}{2m}\cdot X^{\prime\prime}(x)$	$\displaystyle=E_{x}\cdot X(x),$		(27)
	$\displaystyle-\frac{\hbar^{2}}{2m}\cdot Y^{\prime\prime}(y)+mgy\cdot Y(y)$	$\displaystyle=E_{y}\cdot Y(y),$		(28)

where $E_{x}$ and $E_{y}$ are the corresponding energy values in $x-$ and $y-$direction adding up to the total particle energy $E=E_{x}+E_{y}$.

It is obvious that allowed wave functions for the $x-$direction correspond to those of the infinite potential well

$$X_{n}(x)=\sqrt{\frac{2}{L}}\cdot\sin\left(\frac{n\cdot\pi\cdot x}{L}\right),$$

(29)

with $n\in\mathbb{N}^{+}$. The quantization in $x-$direction arises from the boundary conditions $X(0)=0=X(L)$. The associated quantized energy eigenvalues read

$$E_{x,n}=\frac{n^{2}\cdot\pi^{2}\cdot\hbar^{2}}{2mL^{2}}.$$

In order to find solutions to the $y-$direction wave function in Eq. (28) we first define new variables, namely

$$R=\left(\frac{\hbar^{2}}{2m^{2}g}\right)^{\frac{1}{3}},$$

(30)

as well as

$$z=\frac{y}{R}-\epsilon_{y}~{}~{}~{}~{}~{}\mbox{with}~{}~{}~{}~{}~{}\epsilon_{y}=\frac{2mR^{2}\cdot E_{y}}{\hbar^{2}}.$$

(31)

In this new variables Eq. (28) becomes

$$Y^{\prime\prime}(z)-z\cdot Y(z)=0.$$

(32)

The last equation is the well-known Airy-equation. The solutions are the Airy-functions of first $\text{Ai}(z)$ and second kind $\text{Bi}(z)$. General properties and applications of the Airy-functions can e.g. be found in [11]. The general solution for the wave function in $y-$direction is thus given by

$$Y(y)=c_{1}\cdot\text{Ai}\left(\frac{y}{R}-\epsilon_{y}\right)+c_{2}\cdot\text{Bi}\left(\frac{y}{R}-\epsilon_{y}\right).$$

(33)

The boundary conditions at $y=0$ and $y=L$ on $Y(y)$ imply that

	$\displaystyle 0$	$\displaystyle=c_{1}\cdot\text{Ai}\left(-\epsilon_{y}\right)+c_{2}\cdot\text{Bi}\left(-\epsilon_{y}\right),$		(34)
	$\displaystyle 0$	$\displaystyle=c_{1}\cdot\text{Ai}\left(\frac{L}{R}-\epsilon_{y}\right)+c_{2}\cdot\text{Bi}\left(\frac{L}{R}-\epsilon_{y}\right).$		(35)

Since the Airy-functions are transcendental functions there is no closed form solution in general. Because of this, we are now considering two different limiting cases, namely the first where ${\frac{L}{R}\gg\epsilon_{y}}$ and the second, where $\frac{L}{R}-\epsilon_{y}<-1$.

Considering the first case, where $\frac{L}{R}\gg\epsilon_{y}$, simply corresponds to the system, where the energy of the particle in $y-$direction is far beyond the potential energy $mgL$ at the upper boundary. Due to that, one can apply the asymptotic properties of the Airy-function of second kind, namely

$$\lim_{z\rightarrow\infty}\text{Bi}(z)=\infty.$$

(36)

This means for Eq. (35) in order for the total wave function to behave well a the boundary that $c_{2}=0$. The associated boundary condition $Y(0)=0=\text{Ai}(-\epsilon_{y})$ yields the quantized values for $\epsilon_{y}$ and therefore energy eigenvalues in this limit. These energy values can e.g. be determined using a WKB approximation as shown in [9]. The result reads

$$\epsilon_{y,k}=\left(\frac{3\pi}{2}\cdot(k-\frac{1}{4})\right)^{\frac{2}{3}}~{}~{}~{}k\in\mathbb{N}^{+},$$

(37)

and therefore

$$E_{y,k}=\frac{\hbar^{2}}{2mR^{2}}\cdot\left(\frac{3\pi}{2}\cdot(k-\frac{1}{4})\right)^{\frac{2}{3}}.$$

(38)

In order to distinguish the low energy from the high energy regime, we will use the quantum number $k$ for the first, and later on the quantum number $r$, for the other case.

In this limit the appropriate normalisation constant $c_{1}$ of the associated wave function can approximately be determined by the assumption $L\approx\infty$. For the interested reader the integrals involved in the integration are stated in the appendix 7. The total wave function in the low energy limit thus reads

$$Y_{k}(y)=\frac{\text{Ai}\left(\frac{y}{R}-\epsilon_{y,k}\right)}{\sqrt{R\cdot\text{Ai}^{\prime}(-\epsilon_{y,k})^{2}}}.$$

(39)

As a second (and far more interesting case) we consider the situation $\frac{L}{R}-\epsilon_{y,n}<-1$ which is equivalent to the statement that the energy of the particle is larger than the potential energy at the upper boundary. Classically those cases correspond to bounces at $y=0$ and $y=L$. Note that the restriction $<-1$ is chosen in order to keep the errors in the following approximations of Ai and Bi below $1\%$. Solutions to this set of boundary equations (34) & (35) are obtained if

$$\det\begin{pmatrix}\text{Ai}(-\epsilon_{y,n})&\text{Bi}(-\epsilon_{y,n})\\ \text{Ai}\left(\frac{L}{R}-\epsilon_{y,n}\right)&\text{Bi}\left(\frac{L}{R}-\epsilon_{y,n}\right)\end{pmatrix}=0,$$

(40)

is satisfied. To solve the last equation we can approximate the Airy-functions in the considered energy regime via

	$\displaystyle\text{Ai}(x)$	$\displaystyle\approx\frac{\sin\left(\frac{2}{3}(-x)^{\frac{3}{2}}+\frac{\pi}{4}\right)}{\sqrt{\pi}(-x)^{\frac{1}{4}}},$		(41)
	$\displaystyle\text{Bi}(x)$	$\displaystyle\approx\frac{\cos\left(\frac{2}{3}(-x)^{\frac{3}{2}}+\frac{\pi}{4}\right)}{\sqrt{\pi}(-x)^{\frac{1}{4}}}.$		(42)

After some straightforward calculations, including trigonometric identities, the condition (40) can be written compact as

$$\sin\left(\frac{2}{3}\left[\epsilon_{y}^{\frac{3}{2}}-(\epsilon_{y}-\frac{L}{R})^{\frac{3}{2}}\right]\right)=0.$$

(43)

Clearly it is satisfied if

$$\epsilon_{y}^{\frac{3}{2}}-(\epsilon_{y}-\frac{L}{R})^{\frac{3}{2}}=\frac{3r\pi}{2}~{}~{}~{}r\in\mathbb{N}^{+}.$$

(44)

The solutions for the reduced energy eigenvalues $\epsilon_{y}$ depending on $r$ can not be written down in a closed form but due to our restriction $\frac{L}{R}-\epsilon_{y}<-1$ one can Taylor expand the equation up to second order in the variable $\frac{L}{R\epsilon_{y}}=w\approx 0$:

$$\frac{3L}{2R}\cdot\sqrt{\epsilon_{y,r}}-\frac{3L^{2}}{8R^{2}\sqrt{\epsilon_{y,r}}}\approx\frac{3r\pi}{2}.$$

(45)

Solving this equation and inserting the expression for $R$ results in the approximate quantized values for $\epsilon_{y}$

$$\epsilon_{y,r}=\frac{r^{2}\pi^{2}R^{2}}{4L^{2}}\cdot\left[1+\sqrt{1+\frac{L^{3}}{R^{3}\pi^{2}r^{2}}}\right]^{2},$$

(46)

or in terms of the quantized energy eigenvalues

$$E_{y,r}=r^{2}\cdot\frac{E_{1}}{4}\cdot\left[1+\sqrt{1+\frac{mgL}{E_{1}\cdot r^{2}}}\right]^{2}.$$

(47)

Here $E_{1}=\frac{\pi^{2}\hbar^{2}}{2mL^{2}}$ is the ground state energy of a free particle of mass $m$ inside an infinite potential well of length $L$. In the large $r$ limit the asymptotic behavior of the energy equals that of a free particle in the infinite potential well with energy eigenvalues $E_{y,r}\approx E_{1}\cdot r^{2}$. This is the analogous quantum behavior as seen in the classical limit. Note that due to the restriction $\frac{L}{R}-\epsilon_{y,r}<-1$ the allowed values for $r$ depend on the chosen parameters $m$ and $L$ of our system. The complete wave function in $y-$direction thus reads

$$Y_{r}(y)=\frac{\text{Ai}\left(\frac{y}{R}-\epsilon_{y,r}\right)-\frac{\text{Ai}_{-}}{\text{Bi}_{-}}\cdot\text{Bi}\left(\frac{y}{R}-\epsilon_{y,r}\right)}{\sqrt{R\cdot\left[\left(\text{Ai}_{-}^{\prime}-\frac{\text{Ai}_{-}}{\text{Bi}_{-}}\cdot\text{Bi}_{-}^{\prime}\right)^{2}-\left(\text{Ai}_{+}^{\prime}-\frac{\text{Ai}_{-}}{\text{Bi}_{-}}\cdot\text{Bi}_{+}^{\prime}\right)^{2}\right]}}.$$

(48)

The correct normalization constant in terms of derivatives of Ariy-functions has been determined using the integral representations for the Airy-functions displayed in the appendix 7. In order to write the wave function in a compact form we introduced the short hand notation in the last equation, where the index minus refers to the value $-\epsilon_{y,r}$ in the argument of the function and similarly the index plus for $\frac{L}{R}-\epsilon_{y,r}$.

An analogous qualitative plot for the probability density in $y-$direction in the high energy case (approaching the same behavior as a free particle in the infinite potential well) is shown in Figure 7. Note that the fact that the wave function, and therefore the probability density, is not exactly zero at $y=L$ results from the minimal error related to the Taylor expansion for the energy eigenvalues.

Note that both graphs show qualitatively the same behavior as in the classical regime shown in Figure 4. For a three-dimensional representation of the total probability density, consult the graphics displayed in appendix 7.

We close this section with a calculation of the position expectation values and related uncertainties. The $x-$direction quantities can easily be derived (see e.g. [2]) and read in terms of the classical results Eqs. (18) & (20)

	$\displaystyle\langle x\rangle_{qm}$	$\displaystyle=\langle x\rangle_{cl},$		(50)
	$\displaystyle\Delta x_{qm}$	$\displaystyle=\Delta x_{cl}\cdot\sqrt{1-\frac{6}{\pi^{2}n^{2}}},$		(51)

depending on the quantum number $n$.

For the $y-$component quantities we consider the low and high energy limit separately again. The results for the low energy (see Eq. (38)) limit, basically corresponding to a particle without upper boundary, have e.g. been calculated in [2] and read

	$\displaystyle\langle y\rangle_{qm}$	$\displaystyle=\frac{2}{3}\cdot\frac{E_{y,k}}{mg},$		(52)
	$\displaystyle\Delta y_{qm}$	$\displaystyle=\frac{2\sqrt{5}}{15}\cdot\frac{E_{y,k}}{mg}.$		(53)

As in the classical case, we see the identical linear dependence if we associate the quantity $E_{y,k}/mg$ with the theoretical reachable maximal height $h$.

We now turn to the high energy limit where the wave function is given by (48). The expectation value $\langle y\rangle_{qm}$ can be calculated in terms of the variable $z$ introduced in Eq. (31) via

$$\langle y\rangle_{qm}=R\epsilon_{y,r}+R^{2}\int_{-\epsilon_{y,r}}^{\frac{L}{R}-\epsilon_{y,r}}z\cdot|Y_{r}(z)|^{2}dz.$$

(54)

Performing the calculation, including the integral identities stated in the appendix 7 and using the boundary conditions for the wave function, it turns out that many terms appearing in the calculation cancel each other out, resulting in the simple expression for the position expectation value

$$\langle y\rangle_{qm}=\frac{2}{3}\cdot\frac{E_{y,r}}{mg}-\frac{L}{3}\cdot\text{Ji}_{+}.$$

(55)

In the last equation, we defined the new function $\text{Ji}_{+}$ analogous to the classical caseby the expression

$$\frac{\left(\text{Ai}_{+}^{\prime}-\frac{\text{Ai}_{-}}{\text{Bi}_{-}}\cdot\text{Bi}_{+}^{\prime}\right)^{2}}{\left(\text{Ai}_{-}^{\prime}-\frac{\text{Ai}_{-}}{\text{Bi}_{-}}\cdot\text{Bi}_{-}^{\prime}\right)^{2}-\left(\text{Ai}_{+}^{\prime}-\frac{\text{Ai}_{-}}{\text{Bi}_{-}}\cdot\text{Bi}_{+}^{\prime}\right)^{2}},$$

(56)

in order to write the formula for the expectation value in a compact form. Note that $\text{Ji}_{+}$ basically is $R^{2}\cdot Y_{r}^{\prime}\left(\frac{L}{R}-\epsilon_{y,r}\right)^{2}$. As in the classical result (see Eq. (21)), we also see that the expected position in the quantum world is bounded from above by $2E_{y,r}/3mg$ and approaches, for large enough energies, the classical limit.

Let us finally consider the uncertainty in $y-$direction. Due to the translational invariance of the mathematical definition of the uncertainty, it follows that it can be calculated in terms of the variable $z$ as

$$\begin{split}\Delta y_{qm}^{2}=R^{3}\int_{-\epsilon_{y,r}}^{\frac{L}{R}-\epsilon_{y,r}}z^{2}\cdot|Y_{r}(z)|^{2}dz\\ -R^{4}\left(\int_{-\epsilon_{y,r}}^{\frac{L}{R}-\epsilon_{y,r}}z\cdot|Y_{r}(z)|^{2}dz\right)^{2}.\end{split}$$

(57)

Using the integral identities involving Ariy-functions displayed in appendix 7 the calculations fo $\Delta y_{qm}$ can be carried out analytically. Remarkably, due to the boundary conditions, many terms in the lengthy calculation cancel out each other. The quantum uncertainty in terms of $\text{Ji}_{+}$ defined in Eq. (56) for the large energy limit thus takes the form

$$\Delta y_{qm}=\frac{2\sqrt{5}}{15}\cdot\frac{E_{y,r}}{mg}\cdot\sqrt{1+\frac{2mgL}{E_{y,r}}\cdot\text{Ji}_{+}-\frac{7}{2}\cdot\left(\frac{mgL}{E_{y,r}}\cdot\text{Ji}_{+}\right)^{2}}.$$

(58)

Again, as for the expectation value, the similarity to the classical result Eq. (23) is evident.

The classical and quantum study of a particle inside a box subjected to a gravitational field provided interesting results. The conditions on classical periodic orbits, including the interesting geometric phenomena happening at the boundary, have been worked out. A detailed discussion on the probability densities in both classical and quantum regime were derived and considered in specific limiting cases. In particular, the structural similarity of the subsequently obtained formulas for position expectation values, as well as uncertainties, yield a deeper insight into the connection between classical obtained results compared to the quantum ones. For future work, it would be interesting to consider the (quantum) dynamics of systems under the influence of a gravitational field, where either the box is tilted by an angle or the boundary domain is not a square but e.g. a triangle or a parabola (referring to some of the classical study already performed e.g. in [4]).

In this first part of the appendix, some integral identities involving Airy-functions in terms of the variables $z$ defined in Eq. (31), which are used to calculate the normalization constant as well as expectation values, are stated. All of these identities have been obtained using Wolfram Alpha.

$$I_{1}(z):=\int\text{Ai}(z)^{2}=z\cdot\text{Ai}(z)^{2}-\text{Ai}^{\prime}(z)^{2}$$

(59)

$$I_{2}(z):=\int\text{Bi}(z)^{2}=z\cdot\text{Bi}(z)^{2}-\text{Bi}^{\prime}(z)^{2}$$

(60)

$$I_{3}(z):=\int\text{Ai}(z)\cdot\text{Bi}(z)=z\cdot\text{Ai}(z)\cdot\text{Bi}(z)-\text{Ai}^{\prime}(z)\cdot\text{Bi}^{\prime}(z)$$

(61)

$$I_{4}(z):=\int z\cdot\text{Ai}(z)^{2}dz=\frac{1}{6}\left(2z^{2}\text{Ai}(z)^{2}-2z\text{Ai}^{\prime}(z)^{2}+2\text{Ai}(z)\text{Ai}^{\prime}(z)\right)$$

(62)

$$I_{5}(z):=\int z\cdot\text{Ai}(z)\cdot\text{Bi}(z)dz$$

$$=\frac{1}{6}\left(2z^{2}\text{Ai}(z)\text{Bi}(z)+\text{Ai}(z)\text{Bi}^{\prime}(z)+\text{Ai}^{\prime}(z)\text{Bi}(z)-2z\text{Ai}^{\prime}(z)\text{Bi}^{\prime}(z)\right)$$

(63)

$$I_{6}(z):=\int z\cdot\text{Bi}(z)^{2}dz$$

$$=\frac{1}{6}\left(2z^{2}\text{Bi}(z)^{2}-2z\text{Bi}^{\prime}(z)^{2}+2\text{Bi}(z)\text{Bi}^{\prime}(z)\right)$$

(64)

$$I_{7}(z):=\int z^{2}\cdot\text{Ai}(z)^{2}dz$$

$$=\frac{1}{5}\left(\left(z^{3}-1\right)\text{Ai}(z)^{2}-z^{2}\text{Ai}^{\prime}(z)^{2}+2z\text{Ai}(z)\text{Ai}^{\prime}(z)\right)$$

(65)

$$I_{8}(z):=\int z^{2}\cdot\text{Ai}(z)\cdot\text{Bi}(z)dz$$

$$=\frac{1}{5}\left(\text{Ai}(z)\left(\left(z^{3}-1\right)\text{Bi}(z)+z\text{Bi}^{\prime}(z)\right)+z\text{Ai}^{\prime}(z)(\text{Bi}(z)-z\text{Bi}^{\prime}(z))\right)$$

(66)

$$I_{9}(z):=\int z^{2}\cdot\text{Bi}(z)^{2}dz$$

$$=\frac{1}{5}\left(\left(z^{3}-1\right)\text{Bi}(z)^{2}-z^{2}\text{Bi}^{\prime}(z)^{2}+2z\text{Bi}(z)\text{Bi}^{\prime}(z)\right)$$

(67)

Find next some three-dimensional plots for the quantum probability density in the low and high energy limit.