# Finite Volume Effects in Self Coupled Geometries

## I Introduction

In investigating the effects of self coupled boundaries in Monte Carlo lattice QCD simulations there is a need for guidance as to the expected size and direction of boundary-induced energy shifts as well as strategies to remove these effects. Starting with parameterizations of (quenched) lattice potentials within the context of a nonrelativistic two-body Schrödinger-Pauli Hamiltonian with relativistic corrections, a recent approach uses these potentials to produce genuine predictions of charmonium and bottomonium energy levels. These predictions are obtained from numerical solutions, using a purely spatial grid, of the wavefunctions plus perturbation theory. Using this approach, one may begin to understand finite size effects numerically by changing the size of the enclosing box, which is relatively easy to do, unlike the Monte Carlo simulations. Although the general case must still be studied numerically, this article shows that the asymptotic energy shifts from distant boundaries can be determined analytically in terms of the unbounded wavefunctions. Where these wavefunctions are unavailable, it provides functional forms to compare with in potential type simulations, and ultimately, directly to the lattice data once box sizes can confidently be said to be in the asymptotic region.

The systems studied here will be described by the Schrödinger equation in the interior of a volume, the surface of which is determined by a single parameter of a general separable curvilinear coordinate system. It will be assumed that the fields on the boundaries are self-coupled in the sense of being there being either periodic or antiperiodic boundary conditions relating the antipodal points at $`𝐫`$ and $`𝐫`$. We will see that periodic (antiperiodic) boundary conditions require Neumann (Dirichlet) boundary conditions for even parity states and Dirichlet (Neumann) boundary conditions for odd parity states for parity invariant potentials. These boundary conditions are in general not the usual ones applied in cubic geometries, which couple opposite walls rather than antipodal points. However, we will see that the asymptotic shifts when opposite walls are coupled are effectively antipodal. Thus the description here should apply to the situation of Refs., which studies confining solutions (linear plus Coulomb) in a periodic box.

We will proceed to a derivation of the energy shift formulas in Sections II and III. The basic energy shift equations for Dirichlet and Neumann boundary conditions will be derived in Section IV. An asymptotic form of the wavefunction will be presented in Section V, which will then allow us to integrate the differential forms, resulting in simple equations for the energies themselves. It will be seen that the asymptotic energy shifts for Dirichlet and Neumann cases are equal in magnitude but opposite in direction. In Section VI we will make some simplifying observations regarding surface integrals in the energy shift formulas and in Section VII we will explicitly consider the usual lattice situation of cubic geometry with self coupled walls. As examples of the use of these formulas, in Section VIII we will study a one dimensional example, the harmonic oscillator, and a three dimensional example, the hydrogen atom. We will close with some general observations regarding qualitative behaviors of energy shifts in the general case as well as some summary comments.

## II Scrödinger Momentum Tensor Derivation

One way to relate boundary induced energy shifts to properties of the Schrödinger wavefunctions is to use the momentum tensor. The Lagrangian desity for the Schrödinger case may be written

$$=\frac{1}{2m}_i\psi ^{}_i\psi +(VE)\psi ^{}\psi .$$
(1)   
The energy, $`E`$, can be viewed as a Lagrange multiplier associated with the normalization constraint, $`d^3x\psi ^{}\psi =1`$. The Schrödinger equation is recovered from the action defined by Eq.(1) by variation of $`\psi ^{}`$. (In the time dependent case one concludes $`\psi ^{}`$ and $`\psi `$ are related by complex conjugation only after independent variations on $`\psi ^{}`$ and $`\psi `$.) The mass parameter $`m`$ should be understood to represent the reduced mass in applications to two body systems.

The momentum tensor for complex fields may be written as

$$T_{ij}=\frac{\delta }{\delta (_i\psi ^{})}_j\psi ^{}+\frac{\delta }{\delta (_j\psi ^{})}_i\psi ^{}\delta _{ij},$$
(2)   
which is explicitly symmetric in $`i`$ and $`j`$. Using Eq.(1), this gives

$$T_{ij}=\frac{1}{2m}[_i\psi _j\psi ^{}+_j\psi _i\psi ^{}]\frac{\delta _{ij}}{2m}_k(\psi ^{}_k\psi ).$$
(3)   
As a check, we notice that

$$_iT_{ij}=\psi ^{}_jV\psi ,$$
(4)   
which is just a form of Ehrenfest’s theorem when integrated over the volume. The continuity equation

$$\frac{g_i}{t}+_jT_{ji}=0,$$
(5)   
is satisfied in the time dependent case, where the momentum density is given by

$$g_i=\frac{i}{2}(\psi ^{}_i\psi \psi _i\psi ^{}),$$
(6)   
assuring momentum conservation. We will calculate energy shifts by evaluating the pressure exerted by the wavefunction on the surface. For this purpose consider the momentum tensor with locally normal indices:

$$T_{nn}=\frac{1}{2m}_n\psi _n\psi ^{}\frac{1}{2m}_T\psi _T\psi ^{}\frac{1}{2m}\psi ^{}^2\psi .$$
(7)   
Here

$$𝐧\psi _n\psi =\frac{1}{h_n}\frac{\psi }{n},$$
(8)   
is a locally normal derivative ($`𝐧`$ is outwardly directed), $`h_n`$ is a possible scale factor, and $`_T`$ are the tranverse derivatives. In the context of integrations over a closed surface,

$$\frac{1}{2m}\psi ^{}^2\psi +\frac{1}{2m}_T\psi _T\psi ^{}\frac{1}{2m}\psi ^{}_n^2\psi ,$$
(9)   
giving

$$T_{nn}=\frac{1}{2m}|_n\psi |^2\frac{1}{2m}\psi ^{}_n^2\psi .$$
(10)   
Antipodal boundary conditions are given by

$`\psi (𝐫)|_s=\pm \psi (𝐫)|_s,`$ (11)   
$`_n\psi (𝐫)|_s=_n\psi (𝐫)|_s.`$ (12)   
where the top signs for the periodic case and the bottom signs give the antiperiodic case. (One plugs in $`𝐫`$ after the normal derivative is taken in (12).) For a parity invariant potential, one has

$$\psi |_s=0,$$
(13)   
by continuity of the wavefunction at the surface for odd parity states, and

$$_n\psi |_s=0,$$
(14)   
by continuity of the first normal derivative of the wavefunction for even parity states in the periodic case. For antiperiodic boundary conditions, the conditions (13) and (14) apply instead to even and odd parity states, respectively. When the states can not be classified by parity, as when an otherwise good parity system is not centered at the origin, Eqs.(11) and (12) can no longer be simplified.

For the Dirichlet condition Eq.(13), one then has

$$T_{nn}=\frac{1}{2m}|_n\psi |^2,$$
(15)   
and

$$T_{nn}=\frac{1}{2m}\psi ^{}_n^2\psi ,$$
(16)   
for Neumann conditions, Eq.(14).

Finally, we relate the energy shift to an integration over $`T_{nn}`$ via

$$\mathrm{\Delta }E\beta (a_0)=P𝑑V=_{\mathrm{}}^{a_0}𝑑a𝑑sh_nP(a),$$
(17)   
where

$$P(a)=T_{nn}(a).$$
(18)   
Thus, in a differential sense one has

$$\frac{d\beta }{da}=𝑑sh_nT_{nn}(a),$$
(19)   
resulting in

$$\frac{d\beta }{da}=\frac{1}{2m}𝑑sh_n|_n\psi |^2,$$
(20)   
in the Dirichlet case and

$$\frac{d\beta }{da}=\frac{1}{2m}𝑑sh_n\psi ^{}_n^2\psi ,$$
(21)   
in the Neumann. Again, in the general case without good parity, the form from Eq.(10) can not be simplified.

These differential energy shift formulas are actually exact. The wavefunction, $`\psi `$, is a shorthand for $`\psi (\beta (a))|_s`$, denoting the exact wavefunction on the boundary determined at finite $`\beta `$, which is itself determined by the boundary parameter, $`a`$, for a given surface. As they stand, these equations are not yet useful because we have no information on the values of the wavefunction and it’s normal derivatives on $`S`$. A hint on how to proceed is provided by an alternate derivation of these results, which will now be presented.

## III Parameter Variation Derivation

Let us consider the variation of the Schrödinger equation,

$$\frac{1}{2m}^2\psi +V\psi =E\psi ,$$
(22)   
with respect to the parameter (not the coordinate), $`a`$, which determines the self coupling surface, assuming the normalization condition, $`d^3x\psi ^{}\psi =1`$, is satisfied. (Usually the following procedure determines the normalization integral given $`E(a)`$, but here we use it in reverse.) We have

$$\frac{1}{2m}\psi ^{}^2\frac{\psi }{a}+\psi ^{}V\frac{\psi }{a}=\frac{d\beta }{da}\psi ^{}\psi +E\psi ^{}\frac{\psi }{a},$$
(23)   
after variation followed by $`\psi ^{}`$ multiplication. On the other hand, complex conjugation of Eq.(22) followed by multiplication by $`\frac{\psi }{a}`$ results in

$$\frac{1}{2m}\frac{\psi }{a}^2\psi ^{}+\frac{\psi }{a}V\psi ^{}=E\frac{\psi }{a}\psi ^{}.$$
(24)   
The difference produces a perfect differential form, which gives

$$\frac{d\beta }{da}=\frac{1}{2m}𝑑s\left(\frac{\psi }{a}_n\psi ^{}\psi ^{}_n\frac{\psi }{a}\right),$$
(25)   
where $`_n`$ are again local normal derivatives. (The order of the partial derivatives in the last term is immaterial since $`a`$ appears only in $`\psi `$.) We then have

$$\frac{d\beta }{da}=\frac{1}{2m}𝑑s\frac{\psi }{a}_n\psi ^{},$$
(26)   
for Dirichlet boundary conditions, and

$$\frac{d\beta }{da}=\frac{1}{2m}𝑑s\psi ^{}\frac{}{a}_n\psi ,$$
(27)   
for the Neumann case. These results may be reconciled with Eqs.(20) and (21) by the use of expansions in the normal curvilinear coordinate, $`n`$, about the boundary. In the Dirichlet case since the wavefunction vanishes on the surface one must have

$$\psi (na)\frac{\psi }{n}|_s,$$
(28)   
resulting in

$$\frac{\psi }{a}|_s=h_n_n\psi |_s.$$
(29)   
In the Neumann case since the first normal derivative vanishes one has

$$_n\psi (na)\frac{}{n}_n\psi |_s,$$
(30)   
giving

$$\frac{}{a}_n\psi |_s=h_n_n^2\psi |_s.$$
(31)   
This reconciles Eqs.(26) and (27) with (20) and (21). In the general case with no good parity, one must proceed differently. One may expand about the boundary the quantity,

$`\left({\displaystyle \frac{_n\psi (𝐫)}{\psi ^{}(𝐫)}}+{\displaystyle \frac{_n\psi (𝐫)}{\psi ^{}(𝐫)}}\right),`$ (32)   
(for $`|\psi |_s`$ nonvanishing) which also vanishes on S. Then one learns that Eq.(19), with the general $`T_{nn}`$ from (10), and Eq.(25) are equivalent when the antipode terms are grouped together inside the surface integral. One also learns that one may replace

$$_n^2\psi |_s\frac{1}{h_n^2}\frac{^2\psi }{n^2}|_s,$$
(33)   
anywhere in these expressions.

## IV Basic Formulas

The use of the $`a`$ parameter allows one to derive a differential equation for the energy shift, $`\beta `$, and simplifies the Neumann expression. We begin with a change of variables from $`a`$ to $`\beta (a)`$:

$$\frac{\psi }{a}=\frac{d\beta (a)}{da}\frac{\psi }{\beta }.$$
(34)   
In the Dirichlet case one may expand the boundary condition, which determines the energy values, in $`\beta `$,

$$0=\psi |_s\psi (0)|_s+\beta \frac{\psi }{\beta }|_s,$$
(35)   
($`\psi (0)`$ specifies the unbounded wavefunction) allowing the approximate replacement

$$\frac{\psi }{\beta }|_s\frac{1}{\beta }\psi (0)|_s.$$
(36)   
When (29), (34) and (36) are are used in (26), the variables $`\beta `$ and $`a`$ separate and we have

$$^{\beta (a_0)}\frac{d\beta }{\beta ^2}2m^{a_0}\frac{da}{𝑑sh_n^1|\psi (0)|^2}.$$
(37)   
This identifies the anayltic form for $`1/\beta (a_0)`$ up to a constant. Thus, asymptotically one has

$$\beta ^D(a_0)\left(2m^{a_0}\frac{da}{𝑑sh_n^1|\psi (0)|^2}\right)^1.$$
(38)   
The Neumann case proceeds similarly, beginning with

$$\frac{}{a}_n\psi =\frac{d\beta (a)}{da}\frac{}{\beta }_n\psi .$$
(39)   
The boundary condition is again expanded in $`\beta `$,

$$0=_n\psi |_s_n\psi (0)|_s+\beta \frac{}{\beta }_n\psi |_s,$$
(40)   
resulting in

$$\frac{}{\beta }_n\psi |_s\frac{1}{\beta }_n\psi (0)|_s.$$
(41)   
As we will see in the next Section under very general conditions for good parity states,

$$\psi (\beta (a))|_s2\psi (0)|_s.$$
(42)   
When (39), (41) and (42) are used in (27), no integration is needed and we simply obtain

$$\beta ^N(a_0)\frac{1}{m}𝑑s_0\psi ^{}(0)_n\psi (0),$$
(43)   
for the energy shift. Note that we have not assummed that the wavefunctions are separable. For a general confining ellipsoidal surface it can be shown that in the asymptotic limit the appropriate scale factor goes to one. In the following we will also see that it is always possible to choose surfaces for which these factors are effectively unity. Therefore we will set $`h_n=1`$ in the following.

## V Alternate Approach

An alternate approach determines the energy shifts from the form of the asypmtotic wavefunction. Assuming that the potential, $`V`$, is spherically symmetric far from the region where the unbounded wavefunction exists (measured by it’s average radius, $`<r>`$, say), we have the approximate wave equation

$$\left[\frac{d^2}{dr^2}\kappa ^2(r)\right]\psi (\beta (a))=0,$$
(44)   
which is solved by

$`\psi (\beta (a))=\psi (0)|_s(e^{_{r_a}^r\kappa (r^{})𝑑r^{}}\pm e^{_{r_a}^r\kappa (r^{})𝑑r^{}}),`$ (45)   
where $`r_a=r|_s`$ and where the upper sign goes with Neumann boundary conditions and the lower goes with Dirichlet and the second term gives a vanishingly small contribution far from the boundary. In the immediate vincinity of the boundary, the unbounded wavefunction may be characterized as

$$\psi (0)\psi (0)|_se^{_{r_a}^r\kappa (r^{})𝑑r^{}}.$$
(46)   
Thus we conclude in the Dirichlet case,

$$_n\psi (\beta (a))|_s2_n\psi (0)|_s,$$
(47)   
giving the alternate expression,

$$\beta ^D(a_0)\frac{2}{m}_{\mathrm{}}^{a_0}𝑑a𝑑s\left|_n\psi \right|^2.$$
(48)   
Similarly, in the Neumann case one has (42) above as well as

$$_n^2\psi (\beta (a))|_s2_n^2\psi (0)|_s,$$
(49)   
giving the alternate form

$$\beta ^N(a_0)\frac{2}{m}_{\mathrm{}}^{a_0}𝑑a𝑑s\psi ^{}(0)_n^2\psi (0).$$
(50)   
The next Section will show that use of the aymptotic forms of the wavefunctions allows the Dirichlet or Neumann energy shift formulas to be integrated and results in<sup>1</sup>

$$\beta ^D(a_0)\beta ^N(a_0)\frac{1}{m}𝑑s_0\psi ^{}(0)_n\psi (0).$$
(51)   
This equation can also be interpreted as the energy shift in systems in higher or lower dimensions as long as the basic asymptotic wave equation (44) holds in the appropriate variable.

## VI Surface Integrals

In the last two Sections we have derived differently appearing formulas for the energy shifts. These will be reconciled in this Section. In the process we will learn about the surface integrals in these formulas in the general situation when the enclosing surface is not spherical.

The reconcilation of (43) and (50) in the Neumann case is accomplished by use of the asymptotic solutions to Eq.(44). The condition $`r_a<r>`$ in each quantum state (or the appropriate statement in one dimension) and the existence of the asymptotic wave equation, Eq.(44), insures that we can ignore any $`r`$ dependence in the unbounded wavefunction, $`\psi (0)`$, when integrating or taking normal partial derivatives, other than from the exponential in Eq.(46). Writing the general asymptotic form of the unbounded wavefunction (still not assuming separability; $`\mathrm{\Omega }`$ represents the angular variables) as

$$\psi (0)=N(r,\mathrm{\Omega })e^{_0^r𝑑r^{}\kappa (r^{})},$$
(52)   
one can show that either (43) or (50) give the common form,

$$\beta ^N(a_0)\frac{1}{m}𝑑s_0\kappa (r_a)_nr|\psi (0)|^2,$$
(53)   
This form may be further reduced because of the presence of the exponential in the surface integral. This will suppress contributions from regions of the surface not in the immediate vicinity of the points of closest approach to the potential center. Let us call $`r_{}`$ this closest distance. Thus, we may pull the $`r`$ dependent quantities outside of the spatial integral, frozen at their values for $`r=r_{}`$,

$$\beta ^N(a_0)\frac{1}{m}\kappa (r_{})\frac{r_{}}{n}𝑑s_0|\psi (0)|^2.$$
(54)   
In the Dirichlet case, Eq.(48) may similarly be reduced to the negative of the right hand side of Eq.(54) using the asymptotic form, Eq.(52). On the other hand, starting from (38) one may write

$$𝑑s|\psi (0)|^2=e^{2_0^r_{}\kappa (r^{})𝑑r^{}}𝑑s|N(r_a,\mathrm{\Omega })|^2e^{2_r_{}^{r_a}\kappa (r^{})𝑑r^{}}.$$
(55)   
The surface integral in (55) can be treated as a constant in the boundary integral in (38) (verified below), again resulting in the negative of the right hand side of Eq.(54) for $`\beta ^D(a_0)`$.

Eq.(54) is still unnecessarily complicated. Since the near points on the boundary will dominate the surface integral, there are three cases that can occur for the locus of these points: sphere, circle and point. In the sphere case one simply has

$$\beta _{sphere}^N(a_0)\frac{1}{m}\kappa (r_{})𝑑s_0|\psi (0)|^2.$$
(56)   
In the case where the near points are a circle, one can construct a local cylindrical surface for the integration, giving

$`{\displaystyle 𝑑s|N(r_a,\mathrm{\Omega })|^2e^{2_r_{}^{r_a}\kappa (r^{})𝑑r^{}}}`$ $``$ $`r_{}{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z|N(r_{},\mathrm{\Omega }_{})|^2e^{2\kappa _{}(r_ar_{})},`$ (57)   
$``$ $`r_{}{\displaystyle _0^{2\pi }}𝑑\varphi |N(r_{},\mathrm{\Omega }_{})|^2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ze^{\frac{\kappa _{}}{r_{}}z^2},`$ (58)   
$`=`$ $`\sqrt{{\displaystyle \frac{\pi r_{}^3}{\kappa (r_{})}}}{\displaystyle _0^{2\pi }}𝑑\varphi |N(r_{},\mathrm{\Omega }_{})|^2.`$ (59)   
This results in

$$\beta _{circle}^N(a_0)\frac{1}{m}\sqrt{\pi r_{}^3\kappa (r_{})}_0^{2\pi }𝑑\varphi |\psi _{}(0)|^2,$$
(60)   
where $`\psi _{}(0)`$ is the wavefunction along the $`r_{}`$ points. Likewise, if the nearest equidistant points are discrete, one can consider a locally flat surface and evaluate

$`{\displaystyle 𝑑s|N(r_a,\mathrm{\Omega })|^2e^{2_r_{}^{r_a}\kappa (r^{})𝑑r^{}}}`$ $``$ $`2\pi {\displaystyle \underset{i}{}}|N_i(r_{},\mathrm{\Omega }_{})|^2{\displaystyle _0^{\mathrm{}}}𝑑\rho \rho e^{2\kappa _{}(r_ar_{})},`$ (61)   
$``$ $`2\pi {\displaystyle \underset{i}{}}|N_i(r_{},\mathrm{\Omega }_{})|^2{\displaystyle _0^{\mathrm{}}}𝑑\rho \rho e^{\frac{\kappa _{}}{r_{}}\rho ^2},`$ (62)   
$`=`$ $`{\displaystyle \frac{\pi r_{}}{\kappa (r_{})}}{\displaystyle \underset{i}{}}|N_i(r_{},\mathrm{\Omega }_{})|^2,`$ (63)   
where the sum is over the equidistant points. This results in

$$\beta _{points}^N(a_0)\frac{\pi r_{}}{m}\underset{i}{}|\psi _i(0)|^2.$$
(64)   
This last case is the appropriate one for a cubic box enclosure with antipodal boundary conditions.

## VII Boxes with Self Coupled Sides

The antipodal boundary conditions considered here are not the ones used in Monte Carlo lattice simulations since they distinguish one spatial point as the “center”. These simulations generally use periodic or antiperiodic couplings on opposite box sides. Nevertheless, the antipodal conditions are effectively implemented in such simulations because of the local nature of the surface integrals. Let us see how this comes about.

In general one should consider both continuity of $`\psi `$ and $`_n\psi `$ on the surfaces in constructing the energy shift formulas for boxes with self coupled sides. Let us consider the energy shift from the $`x=\pm a`$ surfaces of the box for a parity invariant, centered potential. The asymptotic wavefunction with antipodal boundary conditions in this case can be written

$$\psi (\beta (a))\psi (r_a,+)e^{_{r_a}^r\kappa (r^{})𝑑r^{}}\pm \psi (r_a,)e^{_{r_a}^r\kappa (r^{})𝑑r^{}},$$
(65)   
where the upper sign is for periodic and the lower is for the antiperiodic case, and where

$$\psi (r_a,+)N(r_a,\mathrm{\Omega })e^{_0^{r_a}\kappa (r^{})𝑑r^{}},$$
(66)   
gives the general form of the unbounded wavefunction with $`r`$ replaced by $`r_a=\sqrt{a^2+y^2+z^2}`$. $`\psi (r_a,)`$ is the same as above except with $`\mathrm{\Omega }`$ reflected through the $`x=0`$ plane. The second term on the right of Eq.(65) will be exponentially suppressed in the interior of the box. Now using the general results, Eqs.(10) and (19), and integrating over $`a`$ using the asymptotic form (65), we have

$$\beta ^{box}(a_0)\frac{1}{m}𝑑s_0\kappa (r_a)_nr\mathrm{Re}(\psi ^{}(r_a,+)\psi (r_a,)).$$
(67)   
($`S_0`$ indicates both sides of the box.) We will get the full energy shift by integrating over the other sides of the box in a similar manner. Because of the local nature of these integrations, the contributing portions of the box sides are almost antipodal, and the factors $`\kappa (r_a)`$ and $`_nr`$ can come outside of the surface integral. Then, for wavefunctions with good parity,

$$\psi _{}(r_a,+)=\pm \psi _{}(r_a,),$$
(68)   
Eq.(67) in the periodic case just gives the results of the last Section for Neumann and Dirichlet boundary conditions, respectively. Antiperidicity just reverses these results. Near the box edges the representation (65) fails to be correct, but the local nature of the surface integrations allows it as an approximation. When the potential is not centered or the states do not have good parity, similar considerations show that periodic and antiperiodic energy shifts are negative and equal, at least in the case where one point is much closer to the force center than the other.

## VIII Examples

We will consider two examples of the use of the energy shift formulas, the one dimensional harmonic oscillator and the hydrogen atom.

For the harmonic oscillator, the asymptotic form of the wavefunctions are given by ($`\overline{x}(m\omega )^{1/2}`$; $`\stackrel{~}{x}x_0/\overline{x}`$; $`x_0`$ is the boundary distance from the origin))

$$\psi _N(\stackrel{~}{x})\frac{\stackrel{~}{x}^Ne^{\frac{1}{2}\stackrel{~}{x}^2}}{\sqrt{\overline{x}N!(2)^N\sqrt{\pi }}}.$$
(69)   
Use of (69) in Eq.(51) (the “surface” integral simply supplies a factor of two) gives immediately that

$$\beta _N^D(x_0)\frac{(2)^{N+1}\stackrel{~}{x}^{2N+1}e^{\stackrel{~}{x}^2}}{m\overline{x}^2N!\sqrt{\pi }}.$$
(70)   
On the other hand, the explicit energy shifts for the harmonic oscillator may be found from the general forms of the solutions, proportional to confluent hypergeometic functions, using the boundary condition on the wavefunction. The form of the wavefunction is

$$\psi ^{even}(\stackrel{~}{x})e^{\frac{1}{2}\stackrel{~}{x}^2}F(\frac{\lambda +1}{4}|\frac{1}{2}|\stackrel{~}{x}^2),$$
(71)   
for the even parity states, and the energies are $`E=\lambda /2m\overline{x}^2`$. Expanding $`\lambda `$ as

$$\lambda =(4N+1)+2m\overline{x}^2\mathrm{\Delta }E,$$
(72)   
($`N=0,1,2,\mathrm{}`$) one has in the Dirichlet case , $`\psi ^{even}(\stackrel{~}{x})=0`$, that

$$\mathrm{\Delta }E_N^D(x_0)\frac{F(N|\frac{1}{2}|\stackrel{~}{x}^2)}{2m\overline{x}^2\frac{}{\lambda }F(N\frac{\lambda }{4}|\frac{1}{2}|\stackrel{~}{x}^2)|_{\lambda =0}}.$$
(73)   
Asymptotically,

$$F(N|\frac{1}{2}|\stackrel{~}{x}^2)\frac{(1)^N}{\left(\frac{1}{2}\right)_N}\stackrel{~}{x}^2,$$
(74)   
where

$$(a)_N\underset{i=0}{\overset{N1}{}}(a+i)=\frac{\mathrm{\Gamma }(N+a)}{\mathrm{\Gamma }(a)}.$$
(75)   
The method of taking derivatives of confluent hypergeometric functions in Ref. is used, neglecting polynomials in $`x_0`$ up to order $`2N`$ in $`\frac{}{\lambda }F(N\frac{\lambda }{4},\frac{1}{2},\frac{1}{2}\stackrel{~}{x}^2)|_{\lambda =0}`$. After a shift in the summand and considerable reduction, this can be put into the form

$$\mathrm{\Delta }E_N^D(x_0)\frac{2\stackrel{~}{x}}{N!}\left[\underset{\tau =2N+1}{\overset{\mathrm{}}{}}\frac{\stackrel{~}{x}^{2\tau 4N+1}}{(N+\frac{1}{2})_{\tau 2N}(\tau 2N)_{N+1}}\right]^1.$$
(76)   
In dropping the polynomial contribution, we must have that the major contribution to the infinite sum in (76) occurs for $`\tau 2N`$. Since $`<N|x^2|N>=\overline{x}^2(N+\frac{1}{2})`$ in the unbounded state N and because the major contribution in the infinite sum actually occurs for $`\tau \stackrel{~}{x}^2`$, we have that $`x_0^2\mathrm{\hspace{0.17em}2}<N|x^2|N>`$ for these asymptotic forms to hold, as would be expected.

One notices that the Eq.(70) is odd in $`x_0`$, while (76) is even, so a term by term comparison of the two expressions (expanding $`e^{\stackrel{~}{x}^2}=(e^{\stackrel{~}{x}^2})^1`$ in a power series) is not possible. However, the asymptotic limit of the two sums are identical. This can be established by using

$`(N+{\displaystyle \frac{1}{2}})_{\tau 2N}(\tau 2N)_{N+1}`$ $``$ $`{\displaystyle \frac{(\frac{1}{2})_{\tau +1}2^N}{(2N1)!!}},`$ (77)   
$`{\displaystyle \underset{\tau =2N+1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\stackrel{~}{x}^{2\tau +1}}{(\frac{1}{2})_{\tau +1}}}`$ $``$ $`\sqrt{\pi }e^{\stackrel{~}{x}^2}.`$ (78)   
Eq.(77) holds for $`\tau 2N`$ and Eq.(78) can be justified by using Eq.(75), replacing the sum on $`\tau `$ with an integral, and doing the shift $`\tau \tau \frac{1}{2}`$ (neglecting a polynomial of order $`4N+1`$ in $`\stackrel{~}{x}`$). Finally, using (77) and (78) in (70), we now have that

$$\mathrm{\Delta }E_N^D(x_0)=\beta _{2N}^D(x_0),$$
(79)   
as is appropriate for the even parity states. Going back to the general form of the wavefunction, one may also verify that asymptotically $`\mathrm{\Delta }E_N^N(x_0)=\mathrm{\Delta }E_N^D(x_0)`$ from the Neumann boundary condition, $`_x\psi ^{even}(x)|_{x=x_0}=0`$.

For our second example we consider a hydrogen atom bounded by a sphere. The asymptotic form of the radial wavefunctions for the hydrogen atom are ($`\stackrel{~}{r}r_0/Na_0`$; “$`a_0`$” is the Bohr radius)

$$R_{NL}(\stackrel{~}{r})\left(\frac{2}{Na_0}\right)^{3/2}\frac{(2\stackrel{~}{r})^{N1}e^{\stackrel{~}{r}}}{\sqrt{2N(N+L)!(NL1)!}}.$$
(80)   
This immediately implies the Dirichlet energy shifts

$$\beta _{NL}^D(r_0)\frac{1}{ma_0^2N^3}\frac{(2\stackrel{~}{r})^{2N}e^{2\stackrel{~}{r}}}{(N+L)!(NL1)!},$$
(81)   
from Eq.(51). Obviously, the surface integral was trivial in this case.

The explicit form of the energy shifts for the hydrogen atom have previously been evaluated in Ref. in the Dirichlet case. This was done by expanding the general form of the radial wavefunction

$$R(r)e^{\frac{r}{na_0}}r^LF(Ln+1|2L+2|\frac{2r}{na_0}),$$
(82)   
with energy

$$E=(2ma_0^2n^2)^1,$$
(83)   
to first order in a Taylor series in the energy shift,

$$nN+ma_0^2N^3\mathrm{\Delta }E,$$
(84)   
($`N=1,2,3,\mathrm{}`$) and then solving for $`\mathrm{\Delta }E_{NL}^D`$ with boundary condition $`R(r_0)=0`$. One obtains after reductions,

$$\mathrm{\Delta }E_{NL}^D(r_0)\left[ma_0^2N^3(N+L)!(NL1)!\underset{\tau =2N+1}{\overset{\mathrm{}}{}}\frac{(\tau 2N1)!(2\stackrel{~}{r})^{\tau 2N}}{(\tau LN1)!(\tau +LN)!}\right]^1,$$
(85)   
using the asymptotic form of the wavefunctions and neglecting a finite polynomial series of order $`NL1`$ in $`r_0`$ in the denominator. For $`\tau N`$ or $`L`$ (guaranteed for $`\stackrel{~}{r}N^2`$), one has

$$\frac{(\tau 2N1)!}{(\tau LN1)!(\tau +LN)!}\frac{1}{\tau !},$$
(86)   
which results in

$$\mathrm{\Delta }E_{NL}^D(r_0)=\beta _{NL}^D(r_0),$$
(87)   
when another polynomial of order $`2N`$ in $`r_0`$ is neglected in forming the exponential, $`e^{2\stackrel{~}{r}}`$. Once again, one may verify that the Neumann shift is just the negative of the Dirichlet one.

## IX Comments and Summary

A number of comments may be made about energy shifts in general from the above formulas.

$``$ Averaging over Dirichlet and Neumann boundary conditions by averaging over periodic and antiperiodic boundary conditions in lattice simulations, is an effective way of removing the asymptotic energy shifts for nonrelativistic bound systems with good parity. This point can and should be tested numerically in potential simulations. It is the assumption of good parity which allows this statement to be made, and is not true for antipodal boundary conditions in general.

$``$ The exact formula, Eq.(20), shows that the energy shift is never negative in the case of Dirichlet boundary conditions, which corresponds to the usual meaning of a bounded system. This makes sense from the point of view of the the extra kinetic energy generated from the Heisenberg uncertainty principle.

$``$ The exact differential energy shift in the case of Neumann boundary conditions for a central potential bounded by a sphere can be written as

$$\frac{d\beta }{da}=𝑑s\psi ^{}(\beta (a))\left[V(a)E(a)+\frac{L(L+1)}{2ma^2}\right]\psi (\beta (a)).$$
(88)   
Assuming the right hand side is dominated by the (positive) spherical confinement potential, $`V(a)`$, at large $`a`$, one immediately sees that the asymptotic energy shifts are negative as $`a`$ is decreased. For purely Coulombic systems the right hand side is instead dominated by $`E(a)`$, which is negative, leading again to asymptotic negative energy shifts. At extremely small $`a`$, one expects the positive Neumann box energies, $`E1/a^2`$ to dominate, leading to positive energy shifts. Thus, in the simplest scenario one expects a minimum in the energy at some intermediate value of $`a`$. This seems to be what is seen in the numerical solutions corresponding to Neumann boundary conditions (periodic even parity states) in Ref., although the confining geometry is a cube, not a sphere. Actual proof of this scenario in spherical and other geometries however is difficult and more complicated types of behaviors can not be ruled out.

As illustrations of the derived formulas we considered two examples of centered potentials: a one dimensional harmonic oscillator and a spherically bounded three dimensional hydrogen atom. The asymptotic energy shifts were calculated by the use of the derived formulas and verified from the known forms of the wavefuctions for the even parity oscillator and general quantum number Coulomb systems. Although we have not illustrated the more general case of a confining nonspherical boundary, formulas have been given for the asymptotic surface integrals encountered. The common situation of a confining box with self coupled sides is also amenable to numerical calculation using a formula of Section VII. Thus, the formalism here may readily be applied to lattice simulations of nonrelativistic bound systems in the usual box geometry. When applied to truly periodic systems, not self coupled as in lattice studies, these results actually estimate the location of one edge of the Bloch energy bands which emerge.

## Acknowledgments

This work originated in a conversation with P. Boyle regarding the results in Ref.. The author gratefully acknowledges support from NSF Grant No. PHY-9722073, the Department of Physics and Astronomy, University of Kentucky, as well as from the Special Research Centre for the Subatomic Structure of Matter, University of Adelaide.

FOOTNOTE

1. An alternate unintegrated expression, coming from the Dirichlet expression, Eq.(38), is

$`\beta ^D(a_0){\displaystyle \frac{1}{2m}}{\displaystyle 𝑑s\left[^{a_0}\frac{da}{|\psi (0)|^2}\right]^1}.`$ (89)