The energy eigenfunctions of Schrödinger’s equation for a two-dimensional square-well potential with infinitely high walls are:

φ_{n,p}(x,y,t) |
= | (2/√LM) sin(nπx/L) sin(pπy/M) exp(–2πiE_{n,p}t/h) |
(1) |

where *L* and *M* are the dimensions of the well in the *x*– and *y*-directions (*x* and *y*
being zero at one of the corners of the well), *n* and *p* are integers greater than
or equal to 1, specifying the number of half-wavelengths that fit across the
well in each direction, *h* is Planck’s constant, and *E*_{n,p} is the energy of the
eigenfunction, given by:

E_{n,p} |
= | (h^{2}/8m)((n/L)^{2} + (p/M)^{2}) |
(2) |

where *m* is the mass of the particle. The applet computes the superposition of
these energy eigenfunctions for the current set of amplitudes,
*a*_{n,p}:

φ_{total}(x,y,t) |
= | Σ_{n,p}a_{n,p}φ_{n,p}(x,y,t) |
(3) |

on a grid of *x*, *y* values, and interpolates to find the contours of |φ_{total}|^{2}. (The particle mass and the value for Planck’s constant are chosen
to produce interesting effects on the length and time scales involved.)

To ensure that the probability of the ball lying in each goal is independent of
the grid size, this is computed from first principles (rather than just adding
up the probabilities for the grid points within the goals). The integral of
|φ_{total}|^{2} over an arbitrary rectangular region isn’t too difficult to work
out, and the symmetry of the goals then leads to some simplification, but the
result is too messy to quote in full here. Basically, for any given set of
amplitudes the probability consists of a constant plus a number of cosine waves
at beat frequencies between the eigenfunctions — the constant term coming from
the “beat” between eigenfunctions with themselves (*n*=*k* and *p*=*l* in Equation (4)
below), and with other eigenfunctions with identical energies (other cases where
*E*_{n,p}=*E*_{k,l}). Because of the symmetrical placement of the goals, the
coefficients *I*_{n,k}*J*_{p,l} for the beats between some pairs of eigenfunctions is
precisely zero, hence the unchanging probability for the goals in the game’s
starting configuration.

∫∫_{goal}|φ_{total}|^{2}dx dy |
= | Σ_{n,p}Σ_{k,l}|a_{n,p}a_{k,l}|I_{n,k}J_{p,l}cos(2π(E_{n,p}–E_{k,l})t/h+arg(a_{n,p}*a_{k,l})) |
(4) |

The coefficients *I*_{n,k}*J*_{p,l} in this expression are the components (in
eigenfunction coordinates) of a Hermitian matrix that represents an observable
equal to 1 when the ball is in the goal and 0 when it’s not. The best possible
goal-scoring spectrum (which the applet displays when you shift-click on the
blue spectrum panel) is found by computing the dominant eigenvector of this
matrix (the eigenvector with the greatest eigenvalue). In the
infinite-dimensional case, there’d be an infinite subspace with eigenvalue 1,
but the truncated matrix has a single optimum solution; this maximises the
probability of the ball being in the left-hand goal when *t* is zero (or a
multiple of 100, restoring all the phases) but also yields the same probability
in the right-hand goal when certain phases are reversed (which happens when *t* is
an odd multiple of 50).

The players are modelled as rectangular blips in the potential energy experienced by the ball, too small to have any significant effect on the shape of the wave function when they’re stationary. The matrix for this potential in eigenfunction coordinates is:

H'_{kl,np}(t) |
= | Vδ^{2}/LM [cos((k–n)πx(t)/L) – cos((k+n)πx(t)/L)][cos(( l–p)πy(t)/M) – cos((l+p)πy(t)/M)] |
(5) |

where *V* is the height of the potential, δ is its width, and *x*(*t*) and *y*(*t*) are
the player’s coordinates. When a player moves across the field at a constant
rate, *x*(*t*)=*v*_{x}*t*, *y*(*t*)=*v*_{y}*t*, these matrix elements oscillate; in general, there
are eight harmonic components, all of equal amplitude (give or take a sign),
with frequencies:

f_{kl,np} |
= | (k±n)v_{x}/2L ± (l±p)v_{y}/2M |
(6) |

The time evolution of the eigenfunction amplitudes obeys the set of equations (equivalent to Schiff’s Equation 35.5):

a_{k,l}(t) |
= | (–2π/ih) ∫ Σ_{n,p}H'_{kl,np}(t) a_{n,p}(t) exp(2π(iE_{k,l}–E_{n,p})t/h) dt |
(7) |

The integral in Equation (7) works as a Fourier transform, extracting any
harmonic components of *H*' matching the transition frequency
(*E*_{k,l}–*E*_{n,p})/*h*, or its opposite. So the only significant coupling between
eigenfunctions comes when the value of one of the frequencies in Equation (6)
matches ±(*E*_{k,l}–*E*_{n,p})/*h*.

The applet integrates these equations numerically, including only those cases
where the frequencies match. For the particular values and units used (*h*=1,
*L*=100, *M*=50, *m*=0.00125, *E*_{n,p}=(*n*^{2} + 4*p*^{2})/100), a transition occurs from
φ_{n,p} to φ_{k,l} when:

(k^{2}–n^{2}) + 4(l^{2}–p^{2}) |
= | ±[(k±n)v_{x}/2 ± (l±p)v_{y}] |
(8) |

**Reference:** *Quantum Mechanics* by Leonard I. Schiff, McGraw-Hill, 1968,
sections 9 and 35.