First Spectral Problem: Why Two Decay Conditions Quantize Energy

A second-order ODE has two independent solutions. Requiring decay at one end selects one direction in that two-dimensional solution space. Requiring decay at the other end selects another. A bound state exists when those two directions coincide. This page calculates that condition completely for the harmonic oscillator.

Before starting: you need basic quantum mechanics, differentiation, and elementary trigonometry. The two Airy asymptotic formulas used below are stated explicitly. No previous exact-WKB course is assumed.

First session: allow an estimated 60–90 minutes for the boundary condition, phase matching, exact spectrum, and Exercises 1–2. The numerical laboratory and normalization derivation are extensions. Reader trials have not yet validated this time estimate.

Decay selects a solution at each end

Use dimensionless oscillator units and let $\hbar>0$:

$$-\hbar^2 y''+x^2y=Ey, \qquad y\in L^2(\mathbb R).$$

Here $E$ is the spectral parameter. The coefficient of $y''$ and the potential fix the normalization: this is not the convention $(-\hbar^2 y''+x^2y)/2=Ey$.

Let $y_R$ be the solution decaying as $x\to+\infty$, and $y_L$ the solution decaying as $x\to-\infty$, each with a fixed leading normalization. A nonzero solution decays at both ends precisely when

$$D(E)=\Wr[y_L,y_R]=0, \qquad \Wr[u,v]=uv'-u'v.$$

The Wronskian is constant in $x$ because the equation has no $y'$ term. If the determinant vanishes, the two solutions are proportional; that single solution meets both boundary conditions. Multiplying either endpoint solution by a nonzero constant changes $D$ but preserves its zeros.

For $E>0$, the turning points are $x_\pm=\pm\sqrt E$. Between them, $p(x)=\sqrt{E-x^2}$ is real and solutions oscillate. Outside them, solutions have growing and decaying exponential behavior.

Each turning point supplies a phase

Near the right turning point, the potential is linear to first order: $x^2-E\simeq2\sqrt E(x-x_+)$. After rescaling, the local equation is the Airy equation $Y''(\xi)=\xi Y(\xi)$. Its decaying solution has the asymptotic forms

$$\begin{aligned} \Ai(\xi)&\sim\frac{\xi^{-1/4}}{2\sqrt\pi} \ee^{-2\xi^{3/2}/3},&&\xi\to+\infty,\\ \Ai(-r)&\sim\frac{r^{-1/4}}{\sqrt\pi} \sin\left(\frac23 r^{3/2}+\frac\pi4\right),&&r\to+\infty. \end{aligned}$$

The phase $\pi/4$ belongs to matching the decaying solution through a simple turning point. Applying this connection at both ends gives the leading interior forms, up to nonzero amplitudes,

$$\begin{aligned} y_L(x)&\propto p(x)^{-1/2} \sin\left(\frac1\hbar\int_{x_-}^{x}p(u)\,\dd u+\frac\pi4\right),\\ y_R(x)&\propto p(x)^{-1/2} \sin\left(\frac1\hbar\int_x^{x_+}p(u)\,\dd u+\frac\pi4\right). \end{aligned}$$

The approximation is used away from the turning points and matched to the local Airy solutions near them. Its usual controlled limit keeps the turning points separated as $\hbar$ becomes small. We will verify the resulting oscillator spectrum exactly, including low levels where that limiting argument alone would not suffice.

Match the two oscillatory solutions

Set

$$\theta(x)=\frac1\hbar\int_{x_-}^{x}p(u)\,\dd u, \qquad I(E)=\frac1\hbar\int_{x_-}^{x_+}p(u)\,\dd u.$$

Writing $\alpha=\theta+\pi/4$, the right-hand sine becomes

$$\sin(I-\theta+\pi/4) =\cos I\cos\alpha+\sin I\sin\alpha.$$

The left-decaying solution contains only $\sin\alpha$. The two solutions are proportional throughout the interior when the independent $\cos\alpha$ component disappears:

$$\cos I=0, \qquad I=\pi\left(n+\frac12\right), \quad n=0,1,2,\ldots.$$

The action integral is the area of a semicircle:

$$\int_{-\sqrt E}^{\sqrt E}\sqrt{E-x^2}\,\dd x =\frac{\pi E}{2}.$$

Thus the leading WKB condition predicts $E_n=(2n+1)\hbar$. The half-integer shift has been derived from the two connection phases; it was not inserted into the action integral. For a generic potential, this argument gives a leading approximation and additional corrections can change the spectrum.

An exact calculation for the oscillator

Put

$$z=\sqrt{\frac2\hbar}\,x, \qquad \nu=\frac{E}{2\hbar}-\frac12.$$

The equation becomes $y_{zz}+(\nu+\frac12-z^2/4)y=0$. Define $D_\nu(z)$ as the parabolic-cylinder solution with $D_\nu(z)\sim z^\nu\ee^{-z^2/4}$ for positive large $z$. The endpoint solutions are then

$$y_R(x)=D_\nu\left(\sqrt{2/\hbar}\,x\right), \qquad y_L(x)=D_\nu\left(-\sqrt{2/\hbar}\,x\right).$$

Their exact Wronskian is

$$D(E)=-2\sqrt{\frac\pi\hbar}\, \frac{1}{\Gamma\left(\frac12-\frac{E}{2\hbar}\right)}.$$

The reciprocal gamma function has zeros at $0,-1,-2,\ldots$. Therefore

$$D(E_n)=0 \quad\Longleftrightarrow\quad E_n=(2n+1)\hbar,\qquad n\geq0.$$

This establishes the exact result, independently of the leading WKB argument. Agreement in this solvable model does not make the leading rule exact for an anharmonic potential.

Derive the Wronskian normalization

The parabolic-cylinder values at zero are

$$D_\nu(0)=\frac{2^{\nu/2}\sqrt\pi}{\Gamma((1-\nu)/2)}, \qquad D_\nu'(0)=-\frac{2^{(\nu+1)/2}\sqrt\pi}{\Gamma(-\nu/2)}.$$

The prime here means differentiation with respect to $z$. At $x=0$, the left solution has derivative $-\sqrt{2/\hbar}D_\nu'(0)$ and the right solution has the opposite sign. Hence

$$\Wr_x[y_L,y_R] =2\sqrt{2/\hbar}\,D_\nu(0)D_\nu'(0) =-2\sqrt{\pi/\hbar}\,\Gamma(-\nu)^{-1},$$

where the gamma duplication formula gives the final simplification. As a sign check, $E=0$ gives $D(0)=-2/\sqrt\hbar$. Interchanging $y_L,y_R$ reverses this sign while preserving all zeros.

Compute the spectrum by a different method

Download weber-starter-check.py. It first checks the exact Wronskian at several points away from its zeros. It then independently solves a finite-difference eigenvalue problem on $[-8\sqrt\hbar,8\sqrt\hbar]$, with zero endpoint values. The matrix uses only the differential operator, not its known eigenvalues.

For a grid spacing $\Delta x$, the interior matrix entries are

$$H_{jj}=\frac{2\hbar^2}{\Delta x^2}+x_j^2, \qquad H_{j,j\pm1}=-\frac{\hbar^2}{\Delta x^2}.$$

The supported baseline is Python 3.10 or later with mpmath, NumPy, and SciPy. From the directory containing the downloaded file, run

python3 -m pip install mpmath==1.3.0 numpy==2.0.2 scipy==1.13.1
python3 weber-starter-check.py --json

At $\hbar=1$, the program uses 1200, 2400, and 4800 intervals and extrapolates the leading squared-spacing error. It also increases the endpoint distance at fixed grid spacing. One run gave

Level	Exact energy	Extrapolated grid energy
0	1	1.00000000003
1	3	3.00000000003
2	5	5.00000000007
3	7	7.00000000021

The largest spectral difference was about $2.2\times10^{-10}$, the extrapolation changed by $2.6\times10^{-9}$ between refinements, and the endpoint change was below $7\times10^{-11}$. The declared pass tolerance is $2\times10^{-7}$ in units of $\hbar$; the finer agreement is an observed result, not a rigorous error bound. Last digits can vary with the numerical library and platform.

Try --half-width 2: a small interval can give a stable answer to the wrong boundary problem. The check must fail at the default tolerance. Increasing arithmetic precision alone cannot remove this domain error.

Exercises with progressive help

1. Find the missing phase

Suppose you replace both Airy phases by zero. What matching rule follows, and why does it give the wrong ground-state energy?

Hint 1

Compare $\sin\theta$ with $\sin(I-\theta)$.

Hint 2

The coefficient of $\cos\theta$ must vanish. Compare its zeros with those of $\cos I$ in the correct calculation.

Solution

The incorrect condition is $\sin I=0$, giving the wrong offset $E=2n\hbar$ at positive energies. Formally extending this rule to $n=0$ also permits $E=0$, where the two turning points merge and the preceding simple-turning-point argument no longer applies. Both endpoint-decaying solutions carry a turning-point phase, so the omitted phases change the connection condition. The exact gamma-function Wronskian is nonzero at zero energy.

2. Transfer: a different Hamiltonian normalization

For $\widetilde H=(-\hbar^2\partial_x^2+x^2)/2$, find the energies and explain whether its bound-state functions differ from those above.

Hint

Multiply $\widetilde H y=\widetilde E y$ by two before reusing a formula.

Solution

The ODE has $E=2\widetilde E$, so $\widetilde E_n=(n+1/2)\hbar$. The functions and endpoint conditions are the same. An apparent factor-of-two disagreement can therefore be an operator-convention difference rather than a mathematical error.

3. Transfer: change one boundary condition

Restrict the oscillator to $x\geq0$, require decay at infinity, and impose $y(0)=0$. Which full-line levels remain? What changes for $y'(0)=0$?

Hint 1

The admissible function at infinity is still $D_\nu(\sqrt{2/\hbar}x)$. Use its value or derivative at zero.

Hint 2

For Dirichlet data, the reciprocal gamma argument is $(1-\nu)/2$. For Neumann data, it is $-\nu/2$.

Solution

Dirichlet data select $\nu=1,3,5,\ldots$, the odd full-line states, with $E_k=(4k+3)\hbar$. Neumann data select $\nu=0,2,4,\ldots$, the even states, with $E_k=(4k+1)\hbar$, in both cases $k=0,1,2,\ldots$. The equation alone did not determine the spectrum; the domain and boundary condition changed it.

Continue from the solved problem

Before moving on, explain why two endpoint conditions select discrete energies, derive the half-integer shift without memorizing it, and say which part of the argument is approximate and which is exact.

Continue to formal WKB recursion for higher corrections, exact boundary quantization for the summed connection problem, or the pure quartic oscillator to see a problem where the harmonic-oscillator coincidence no longer suffices. The quartic calculation introduces nonlinear integral equations; use Route B’s preparation and laboratory sessions to approach it in stages.

The source formulas are the Airy asymptotics, the parabolic-cylinder equation and values at zero, and the gamma duplication formula.