# Finding expansion coefficient of a 3-d Gaussian wave packet

• I
I'm having trouble with trying to find the expansion coefficients of a superposition of a Gaussian wave packet.

First I'm decomposing a Gaussian wave packet
$$\psi(\textbf{r},0) = \frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{(\textbf{r} - \textbf{r}_0)^2}{4\sigma^2} + i\textbf{k}_0 \cdot \textbf{r}\right]$$
into a linear superposition of eigenfunctions of the hydrogen atom, where the initial position and momentum is chosen as ##\textbf{r}_0## and ##\textbf{p}_0 = \hbar \textbf{k}_0##, and also a spatial width ##\sigma## such that the initial wave packet ##\psi(\textbf{r},0)## can to a sufficient approximation be superimposed by bound-state eigenfunctions ##\varphi_\textit{nlm}(\textbf{r})##,
$$\psi(\textbf{r},0) = \sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l} b_{\textit{nlm}}\varphi_\textit{nlm}(\textbf{r})$$
Having solved the Schrodinger Equation of a hydrogen atom with a Hamiltonian containing the Coulomb potential, the eigenfunctions were obtained to be a product of the Radial Wave solution and Spherical Harmonics
$$\varphi_\textit{nlm}(\textbf{r}) = R_{nl}(r) Y_{lm}(\theta, \phi)$$
where
$$R_{nl}(r) = \frac{1}{a^{3/2}}\frac{2}{n^2}\sqrt{\frac{(n-l-1)!}{(n+1)!}}\left(\frac{2r}{na}\right)^l e^{-\frac{r}{na}} L^{2l+1}_{n-l-1}\left(\frac{2r}{na}\right)$$
where ##L^{2l+1}_{n-l-1}\left(\frac{2r}{na}\right)## is the Laguerre Polynomial which has the form of
$$L^{k}_{p}\left(x\right) = \sum_{k}^{p} (-1)^s \binom {p+k}{p-s} \frac{x^s}{s!}$$
and
$$Y_{lm}(\theta, \phi) = (-1)^m \sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}} P^{m}_{l} (\cos\theta) e^{im\phi}$$
with the Legendre Polynomial contained within the Associated Legendre Functions
$$P^{m}_{l} (x) = (1 - x^2)^{\frac{m}{2}} \frac{d^m}{dx^m}(P_{l}(x))$$
$$P_{l}(x) = \frac{1}{2^l l! }\frac{d^l}{dx^l}[(x^2 - 1)^l]$$

What I ultimately would like to find is the closed form expansion coefficients ##b_{nlm}## of the eigenfunction superposition, now from my knowledge of finding the expansion coefficients of a single-indexed superposition of a wavefunction, the same method can be applied on this superposition which involves 3 different indices, such that
$$\begin{equation*} \begin{split} \langle\varphi_\textit{nlm}(\textbf{r})|\psi(\textbf{r},0)\rangle &= \langle\varphi_\textit{nlm}(\textbf{r})|\sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l} b_{\textit{nlm}}\varphi_\textit{nlm}(\textbf{r})\rangle\\ &= b_{\textit{nlm}}\sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l}\langle\varphi_\textit{n'l'm'}(\textbf{r})| \varphi_\textit{nlm}(\textbf{r})\rangle\\ &= b_{\textit{nlm}}\sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l}\delta_{n'n}\delta_{l'l}\delta_{m'm}\\ &= b_{nlm} \end{split} \end{equation*}$$

My attempt:
To avoid doing the scalar product directly and having to compute the whole integral, I was thinking of using the decomposition of stationary harmonic plane wave into partial waves, as such
$$e^{\textbf{k}\cdot\textbf{r}} = e^{ikz} = e^{ikr\cos\theta} = \sum^{l=0}_{\infty}(2l+1)i^l j_l(kr) P_l(\cos\theta)$$
where ##i## is the imaginary number, ##\textbf{k}## and ##\textbf{r}## are the wave and position vectors, ##j_l## are the spherical Bessel functions, ##P_l## are Legendre Polynomials. But according to Wikipedia (https://en.wikipedia.org/wiki/Plane_wave_expansion), it can also be rewritten as
$$e^{\textbf{k}\cdot\textbf{r}} = 4\pi \sum_{l=0}^{\infty}\sum_{m=-l}^{l} i^l j_l(kr)Y^m_l(\hat{\textbf{k}})Y^{m*}_l(\hat{\textbf{r}})$$

In theory, this should now mean that I can rewrite the scalar product in terms of these sum of spherical waves, so that

$$\begin{equation*} \begin{split} b_{nlm} &= \langle\varphi_\textit{nlm}(\textbf{r})|\psi(\textbf{r},0)\rangle\\ &= \langle R_{nl}(r) Y_{lm}(\theta, \phi)|\frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{(\textbf{r} - \textbf{r}_0)^2}{4\sigma^2}\right]\cdot 4\pi\sum_{l=0}^{\infty}\sum_{m=-l}^{l} i^l j_l(kr)Y^m_l(\hat{\textbf{k}})Y^{m*}_l(\hat{\textbf{r}})\rangle\\ \end{split} \end{equation*}$$
and so this is where I'm at. I'm not sure how to go from here. I hope I've clarified my question well.

hilbert2
Gold Member
Is this related to "Kepler wavepackets", where an electron wavepacket orbits a nucleus like in classical central-force motion?

Is this related to "Kepler wavepackets", where an electron wavepacket orbits a nucleus like in classical central-force motion?
Yes how did you know?Ultimately I'm trying to investigate the correspondence principle of the hydrogen atom, so in the limit of high quantum number how it'd would collapse from a quantum mechanical regime to a classical one.

hilbert2
Gold Member
Someone I worked with many years ago made his master's thesis about that subject.

There seems to be some related material that can be found with a Google search, but I'm not sure if it's useful to you. I suggest you represent the Gaussian wavepacket in spherical polar coordinates, as the hydrogenic eigenstates are, and then try to calculate the expansion coefficients by integration with Wolfram Mathematica. If that application can't solve the integrals exactly, it's very unlikely that they could be found with pen and paper, either.

Someone I worked with many years ago made his master's thesis about that subject.

There seems to be some related material that can be found with a Google search, but I'm not sure if it's useful to you. I suggest you represent the Gaussian wavepacket in spherical polar coordinates, as the hydrogenic eigenstates are, and then try to calculate the expansion coefficients by integration with Wolfram Mathematica. If that application can't solve the integrals exactly, it's very unlikely that they could be found with pen and paper, either.
My supervisor did his master's thesis on it too but he did it on a Quantum Harmonic Oscillator. Most of the stuff can be found from the book "The Picture Book of Quantum Mechanics" but they didn't show how to find the expansion coefficients. This project is waaay above my knowledge of QM which is why I'm always stuck at every step of progress that I'm trying to make lol.

hilbert2
Gold Member
Based on what I can find about it problem, it doesn't seem to be possible to find a closed form solution for the overlap integral of a gaussian wavepacket and the hydrogenic orbitals. It may be possible for a wavepacket in a 3d harmonic oscillator potential ##V(x,y,z)=\frac{1}{2}k(x^2 + y^2 + z^2 )##.

Based on what I can find about it problem, it doesn't seem to be possible to find a closed form solution for the overlap integral of a gaussian wavepacket and the hydrogenic orbitals. It may be possible for a wavepacket in a 3d harmonic oscillator potential ##V(x,y,z)=\frac{1}{2}k(x^2 + y^2 + z^2 )##.
It seems that the Legendre functions are related to the Spherical Harmonics by $$P_l(\cos\theta) = \sqrt{\frac{4\pi}{2l+1}} Y_{l,0}(\theta, \phi)$$ when m = 0, so the decomposition of stationary harmonic plane wave into partial waves can now be expressed as $$e^{\textbf{k}\cdot\textbf{r}} = e^{ikz} = e^{ikr\cos\theta} = \sum^{\infty}_{l=0}(2l+1)i^l j_l(kr) P_l(\cos\theta) = \sqrt{\frac{4\pi}{2l+1}}\sum^{\infty}_{l=0}(2l+1)i^l j_l(kr) Y_{l,0}(\theta, \phi)$$
which I can then incorporate into my Gaussian wave packet, but how can I go about changing the ##(\textbf{r} - \textbf{r}_0)^2## into spherical polar coordinates? I'm only trying this approach because my supervisor said to "expand wavepackets into plane waves, then expand plane waves into spherical harmonics", and I think hopefully the spherical harmonics would cancel out nicely due to them being orthonormal when I do the scalar product.

Which also sounds right because partial waves decomposition works well with a central potential, which is exactly the case with this problem because the potential in the Hamiltonian is central in the Hamiltonian of the SE that I had to solve earlier.

hilbert2
Gold Member
but how can I go about changing the ##(\textbf{r} - \textbf{r}_0)^2## into spherical polar coordinates?

That's not a problem. If the position vectors are ##\mathbf{r}=(x,y,z)## and ##\mathbf{r}_0=(x_0 ,y_0 ,z_0 )##, then

##(\mathbf{r}-\mathbf{r}_0 )^2 = (x-x_0 )^2 + (y-y_0 )^2 + (z-z_0 )^2##,

and then you just use

##x = r\sin\theta \cos\phi##
##y = r\sin\theta \sin\phi##
##z = r\cos\theta##.

Only the ##x,y,z## are changed to these polar coordinates, the ##x_0 ,y_0 ,z_0## stay as they are, because they're just numbers and not variables.

Dazzabaijan
That's not a problem. If the position vectors are ##\mathbf{r}=(x,y,z)## and ##\mathbf{r}_0=(x_0 ,y_0 ,z_0 )##, then

##(\mathbf{r}-\mathbf{r}_0 )^2 = (x-x_0 )^2 + (y-y_0 )^2 + (z-z_0 )^2##,

and then you just use

##x = r\sin\theta \cos\phi##
##y = r\sin\theta \sin\phi##
##z = r\cos\theta##.

Only the ##x,y,z## are changed to these polar coordinates, the ##x_0 ,y_0 ,z_0## stay as they are, because they're just numbers and not variables.
Just in case that you're still interested in helping me, all of the equations above were followed closely from Chapter 10.1, 10.9, and 10.10 of the book "The Picture Book of Quantum Mechanics 4th edition", so it might be useful if you have a look. Meanwhile I'm trying all the approach I could lol

hilbert2
Gold Member
I found the 3rd edition of that book as a pdf document, and it looks quite interesting. The section about Kepler wavepackets doesn't seem to go very far in explaining how to actually do the calculations.

If you have Wolfram Mathematica available, you can try to calculate the integral of the wavepacket multiplied by a hydrogenic wavefunction ##\psi_{n,l,m_l}(r,\phi ,\theta )## over all 3d space, representing all position variables in spherical polar coordinates. It's likely that it's not possible, and the integrations have to be done by some numerical method.

Just used Mathematica for the first time and had to figure out how to do integration. So doing the whole integral of the scalar product gave me this, not sure if it's right or what to make of it. Any opinions would be appreciated!

#### Attachments

• Integral.png
19.2 KB · Views: 140
That's not a problem. If the position vectors are ##\mathbf{r}=(x,y,z)## and ##\mathbf{r}_0=(x_0 ,y_0 ,z_0 )##, then

##(\mathbf{r}-\mathbf{r}_0 )^2 = (x-x_0 )^2 + (y-y_0 )^2 + (z-z_0 )^2##,

and then you just use

##x = r\sin\theta \cos\phi##
##y = r\sin\theta \sin\phi##
##z = r\cos\theta##.

Only the ##x,y,z## are changed to these polar coordinates, the ##x_0 ,y_0 ,z_0## stay as they are, because they're just numbers and not variables.
So following your advice of converting my wavepacket wavefunction into spherical polar coordinates, I get
$$\begin{equation*} \begin{split} \psi(\textbf{r},0) &= \frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{(\textbf{r} - \textbf{r}_0)^2}{4\sigma^2} + i\textbf{k}_0 \cdot \textbf{r}\right]\\ &= \frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{1}{4\sigma^2}[(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2]\right] \text{exp}\left[i\textbf{k}_0 \cdot \textbf{r}\right]\\ &= \frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{1}{4\sigma^2}[(r\sin\theta\cos\phi - x_0)^2 + (r\sin\theta\sin\phi - y_0)^2 + (r\cos\theta - z_0)^2]\right] \sqrt{\frac{4\pi}{2l+1}}\sum^{\infty}_{l=0}(2l+1)i^l j_l(kr) Y_{l,0}(\theta, \phi)\\ &= \eta(r, \theta, \phi) \sqrt{4\pi}\sum^{\infty}_{l=0}\sqrt{(2l+1)}i^l j_l(kr) Y_{l,0}(\theta, \phi) \end{split} \end{equation*}$$
So then I thought I'd use this expression to do the inner product instead, so I got to

$$\begin{equation*} \begin{split} b_{nlm} &= \langle\varphi_\textit{nlm}(\textbf{r})|\psi(\textbf{r},0)\rangle\\ &= \langle R_{nl}(r) Y_{lm}(\theta, \phi)|\eta(r, \theta, \phi) \sqrt{4\pi}\sum^{\infty}_{l=0}\sqrt{(2l+1)}i^l j_l(kr) Y_{l,0}(\theta, \phi)\rangle\\ \end{split} \end{equation*}$$
Then I realised here it's starting to look like that I should be able to use the orthonormalization relations of Spherical Harmonics, where
$$\int_{\theta = 0}^{\pi}\int_{\phi = 0}^{2\pi} Y^m_l {Y^{m'}_{ l'}}^{*} \sin\theta d\phi d\theta = \delta_{ll'}\delta_{mm'}$$
since ##m## for one of the Spherical Harmonics is already 0, we can get an expression for using that relation for ##l=l'## and ##m=m'=0##, but I can't really do that with ##\eta(r, \theta, \phi)## in place, since I've converted it into using spherical polar coordinates. The way I see it could only work if I can somehow change ##\frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{(\textbf{r} - \textbf{r}_0)^2}{4\sigma^2}\right]## to a scalar function that only depends on ##r##.