Spherical HO basis

The standard spherical HO wave functions, are given by

$$\phi_{Nljm}(r,\theta,\phi,\sigma) \!\! = b^{3/2} F_{Nl}(br)e^{-{\textstyle{\frac{1}{2}}}(br)^2} \sum_{m_lm... ...c{1}{2}}}m_s}Y_{lm_l}(\theta,\phi)\chi_{{\textstyle{\frac{1}{2}}}m_s}(\sigma) ,$$ (78)

where $b$ is the oscillator constant,

$$b=\sqrt{\frac{m\omega}{\hbar}},$$ (79)

and $F_{Nl}(br)=$ are proportional to the standard Laguerre polynomials [7].

To calculate the secondary densities (76), one has to act on the space part of the wave function (79) with the derivative operators $D_{nLM_L}$. This leads to defining the polynomials $F^{km_k}_{nLM_LNlm_l}(br)$ such that

$$D_{nLM_L} b^{3/2} F_{Nl}(br)e^{-{\textstyle{\frac{1}{2}}}(br)^2}Y_{lm_l}(\theta,\phi)$$

$$=b^{n+3/2} \sum_{km_k} F^{km_k}_{nLM_LNlm_l}(br) e^{-{\textstyle{\frac{1}{2}}}(br)^2} Y_{km_k}(\theta,\phi) ,$$ (80)

where $m_k=M_L+m_l$. From the Wigner-Eckart theorem, these polynomials must have the form:

$$F^{km_k}_{nLM_LNlm_l}(br) = {\textstyle{\frac{1}{\sqrt{2k+1}}}} C^{km_k}_{lm_lLM_L} F^{k}_{nLNl}(br),$$ (81)

and our goal is to determine the set of reduced polynomials $F^{k}_{nLNl}(br)$, in terms of which the derivatives of the spherical HO wave functions read

$$D_{nLM_L}\phi_{Nljm}(r,\theta,\phi,\sigma) \!\! = b^{n+3/2}\sum_{km_k}\sum_{m_lm_s} {\textstyle{\frac{1}{\sqrt{2k+1}}}} C^{km_k}_{lm_lLM_L}C^{jm}_{lm_l{\textstyle{\frac{1}{2}}}m_s}$$

$$\times F^{k}_{nLNl}(br)e^{-{\textstyle{\frac{1}{2}}}(br)^2}Y_{km_k}(\theta,\phi)\chi_{{\textstyle{\frac{1}{2}}}m_s}(\sigma) ,$$

$$= -b^{n+3/2}\sqrt{2j+1}\sum_{km_k}\sum_{m_s} \sum_{im_i} C^{im_i}_{LM_Ljm}C^{im_i}_{{\textstyle{\frac{1}{2}}}m_skm_k}$$

$$\times \left\{\begin{array}{rrr} l & L & k \\ i & {\textstyle{\frac{1... ...}{2}}}(br)^2}Y_{km_k}(\theta,\phi)\chi_{{\textstyle{\frac{1}{2}}}m_s}(\sigma) .$$

Explicit form of $F^{k}_{nLNl}(br)$ can be calculated by using the Wigner-Eckart theorem again, namely, Eq. (81) must have the form

$$D_{nLM_L} \phi_{Nlm_l} = \sum_{N'km_k} {\textstyle{\frac{1}{\sqrt{2k+1}}}}C^{km_k}_{lm_lLM_L} \langle\phi_{N'k}\vert\vert D_{nL}\vert\vert\phi_{Nl}\rangle \phi_{N'km_k},$$ (82)

where $\phi_{Nlm_l}$ is the space part of the wave function (79). Then we have polynomials $F^{k}_{nLNl}(br)$ expressed through the Laguerre polynomials as:

$$b^n F^{k}_{nLNl}(br) = \sum_{N'} \langle\phi_{N'k}\vert\vert D_{nL}\vert\vert\phi_{Nl}\rangle F_{N'k}(br),$$ (83)

where the reduced matrix element can be calculated by considering only one matrix element, namely,

$$\langle\phi_{N'k,m_k=0}\vert D_{nL,M_L=0}\vert\phi_{Nl,m_l=0}\rangle = {\textstyle{\frac{1}{\sqrt{2k+1}}}}C^{k0}_{l0L0}\langle\phi_{N'k}\vert\vert D_{nL}\vert\vert\phi_{Nl}\rangle .$$ (84)

Note that the Clebsh-Gordan coefficient $C^{k0}_{l0L0}$ is not zero for angular momenta restricted by the parity conservation, $(-1)^{L+l-k}=1$.

The parity conservation induces specific conditions on the polynomials $F^{k}_{nLNl}(br)$. Indeed, by comparing parities of both sides of Eq. (81) we see that

$$F^{k}_{nLNl}(br) = 0 ,\quad\mbox{for}\quad (-1)^{n+l+k}=-1.$$ (85)

Equivalently, since for all derivative operators we have $(-1)^{n+L}=1$, we see that

$$F^{k}_{nLNl}(br) = 0 ,\quad\mbox{for}\quad (-1)^{L+l+k}=-1.$$ (86)

Since polynomials $F_{Nl}(br)$ are real, phases of polynomials $F^{k}_{nLNl}(br)$ are fixed by those of the derivative operators (31) and spherical harmonics [6],

$$Y_{JM}^*(\theta,\phi) = (-1)^{-M}Y_{J,-M}(\theta,\phi) ,$$ (87)

that is,

$$F^{k*}_{nLNl}(br) = (-1)^{l-k}F^{k}_{nLNl}(br) = (-1)^{n}F^{k}_{nLNl}(br) = (-1)^{L}F^{k}_{nLNl}(br),$$

where the last two equivalent forms result from Eqs (87) and (88).

The $F_{nLNl}^k$ polynomials are calculated using formulas for spherical derivatives (see [6] ) combined with recursion relations for derivatives of Laguerre polynomials. In this way spherical derivatives of one of the basis functions $\phi_{Nljm}\left(\vec{r},\sigma\right) = g_{nl}\left(r\right) C_{lm_l,\frac{1}{2}\sigma}^{jm} Y_{lm}\left( \theta,\phi \right)$ can be expressed as a sum of functions

$$\begin{eqnarray*} \nabla_{1\mu_{1}}..\nabla_{1\mu_{N}}g_{nl}\left(r\right)Y_{lm}... ...}\left(r\right)\right]Y_{l+i_{1}+..+i_{n},m+\mu_{1}+..+\mu_{N}}, \end{eqnarray*}$$

where the $a$ and $b$ coefficients needed for the different orders were derived using symbolic programming. In this way we obtain analytical expressions for all derivatives which can then be calculated with good accuracy.