Fields

Now we are in a position to perform the variation of densities over the wave functions. To this end, we assume that the non-local density matrices in Eq. (40) have the general form of

$$\rho\left(\vec{r}\sigma,\vec{r}\,'\sigma'\right) =\sum_i \phi_i(\vec{r}\sigma)\psi_i(\vec{r}\,'\sigma').$$ (46)

This form allows us to carry out the derivation for several important specific cases simultaneously. Namely, for the standard HF case one has:

$$\psi_i(\vec{r}\,'\sigma') = \phi_i^*(\vec{r}\,'\sigma'),$$ (47)

and the sum runs over the occupied states only, $i=1,\ldots,A$. Similarly, for the BCS case, or for the HFB case in the canonical basis, one has:

$$\psi_i(\vec{r}\,'\sigma') = v_i^2\phi_i^*(\vec{r}\,'\sigma'),$$ (48)

where $v_i^2$ are the occupation factors and the sum runs over the pairing window. For transition densities pertaining to the symmetry restoration, one has:

$$\psi_i(\vec{r}\,'\sigma') = \sum_j O_{ij}^{-1}(\alpha)\phi_j^*(\vec{r}\,'\sigma',\alpha),$$ (49)

Where $\phi_j(\vec{r}\,'\sigma',\alpha)=R(\alpha)\phi_j(\vec{r}\,'\sigma')$ are the wave functions transformed by the symmetry operator $R(\alpha)$ and $O_{ij}(\alpha)$ is the overlap matrix:

$$O_{ij}(\alpha) = \int {\rm d}^3\vec{r}\sum_\sigma \phi_i^*(\vec{r}\sigma,\alpha) \phi_j(\vec{r}\sigma).$$ (50)

Finally, for the RPA amplitudes given by non-hermitian matrix $\tilde\rho_{ij}$ one has

$$\psi_i(\vec{r}\,'\sigma') = \sum_j \tilde\rho_{ij}^*\phi_j^*(\vec{r}\,'\sigma').$$ (51)

To derive the fields, first we recall that the variation of the non-local density, over the wave function reads

$$\delta\rho_{v\mu}\left(\vec{r},\vec{r}\,'\right) = \sum_i\sum_{\sigma'}\frac{\partial\rho_{v\mu}\left(\vec{r},\vec{r}\,'\right)} {\partial\psi_i(\vec{r}\,'\sigma')}\delta\psi_i(\vec{r}\,'\sigma')$$

$$= \sum_i\sum_{\sigma\sigma'} \left\langle \sigma'\left\vert\sigma_{... ...t\vert\sigma\right\rangle \phi_i(\vec{r}\sigma)\delta\psi_i(\vec{r}\,'\sigma').$$ (52)

Therefore, variation of the primary density is given by

$$\delta\rho_{n'L'v'J'M'} = \sum_i\sum_{\sigma\sigma'} \left\{ [K_{n'L'}\sigma_{v';\sigma'\si... ...phi_i(\vec{r}\sigma)\delta\psi_i(\vec{r}\,'\sigma')\right\}_{\vec{r}'=\vec{r}}.$$ (53)

At this point, operators $K_{n'L'}$ mix derivatives acting on the variables $\vec{r}$ and $\vec{r}\,'$, cf. Eq. (18). By using the Wigner-Eckart theorem, we can express them as sums of products of derivatives $D_{mIM_I}$ and $D'_{m'I'M'_I}$, which act on $\vec{r}$ and $\vec{r}\,'$, respectively, that is,

$$K_{n'L'M'_L} = \sum_{mIm'I'} K^{n'L'}_{mIm'I'}\sum_{M_IM_I'} C^{L'M'_L}_{IM_II'M'_I}D_{mIM_I}D'_{m'I'M'_I},$$ (54)

where the order of derivative is conserved,

$$n'=m+m' ,$$ (55)

and the triangle rule of angular momentum coupling must be obeyed. Numerical coefficients $K^{n'L'}_{mIm'I'}$ can be calculated by using methods of symbolic programming. At N$^3$LO, only 91 coefficients $K^{n'L'}_{mIm'I'}$ are needed, so they can easily be precalculated and stored.

We now can insert expressions (54) and (55) into Eq. (45) and remove condition ${\vec{r}'=\vec{r}}$ by adding the integral over $\vec{r}\,'$ and the function $\delta(\vec{r}-\vec{r}\,')$, that is,

$$\delta{\cal E} = \sum_{n'L'v'J'}\int {\rm d}^3\vec{r}\int {\rm d}^3\vec{r}\,'\delta(\vec{r}-\vec{r}\,')$$

$$\times \sum_i\sum_\sigma [\tilde{U}_{n'L'v'J'}(\vec{r}) [K_{n'L'}\sigma_{v';\sigma'\sigma}]_{J'}]_0 \phi_i(\vec{r}\sigma)\delta\psi_i(\vec{r}\,'\sigma').$$ (56)

This allows us to integrate by parts over $\vec{r}\,'$ and transfer the action of $D'_{m'I'M'_I}$ onto the delta function. With all the angular momenta couplings shown explicitly, this gives

$$\delta{\cal E} = \sum_{n'L'v'J'}\int {\rm d}^3\vec{r}\int {\rm d}^3\vec{r}\,'\sum_{mIm'I'}\sum_{M_IM_I'} \left\{{D'}^T_{m'I'M'_I}\delta(\vec{r}-\vec{r}\,')\right\}$$

$$\times \sum_i\sum_{\sigma\sigma'} \sum_{M'_L\mu'}\sum_{M'}C^{J'M'}_{L'M'... ...\textstyle{\frac{(-1)^{J'-M'}}{\sqrt{2J'+1}}}}\tilde{U}_{n'L'v'J',-M'}(\vec{r})$$

$$\times K^{n'L'}_{mIm'I'} C^{L'M'_L}_{IM_II'M'_I}D_{mIM_I}\sigma_{v'\mu';\sigma'\sigma} \phi_i(\vec{r}\sigma)\delta\psi_i(\vec{r}\,'\sigma').$$ (57)

The action of ${D'}^T_{m'I'M'_I}$ onto the delta function can be replaced by that of $(-1)^{m'}D^T_{m'I'M'_I}$, and then the integration by parts over $\vec{r}$ gives,

$$\delta{\cal E} = \sum_{n'L'v'J'}\int {\rm d}^3\vec{r}\int {\rm d}^3\vec{r}\,'\sum_{mIm'I'}\sum_{M_IM_I'} \delta(\vec{r}-\vec{r}\,')$$

$$\times \sum_i\sum_{\sigma\sigma'} \sum_{M'_L\mu'}\sum_{M'}C^{J'M'}_{L'M'... ...^{J'-M'}}{\sqrt{2J'+1}}}}(-1)^{m'}D_{m'I'M'_I}\tilde{U}_{n'L'v'J',-M'}(\vec{r})$$

$$\times K^{n'L'}_{mIm'I'} C^{L'M'_L}_{IM_II'M'_I}D_{mIM_I}\sigma_{v'\mu';\sigma'\sigma} \phi_i(\vec{r}\sigma)\delta\psi_i(\vec{r}\,'\sigma').$$ (58)

In the resulting local integral we request that the operator acting on $\phi_i(\vec{r}\sigma)$ is equal to the mean-field operator, which gives

$$h(\rho) = \sum_{n'L'v'J'}\sum_{mIm'I'}\sum_{M_IM_I'} \sum_{M'_L\mu'} \sum_{M'}C^{J'M'}_{L'M'_Lv'\mu'} {\textstyle{\frac{(-1)^{J'-M'}}{\sqrt{2J'+1}}}}$$

$$\times (-1)^{m'}K^{n'L'}_{mIm'I'} C^{L'M'_L}_{IM_II'M'_I} D_{m'I'M'_I}\tilde{U}_{n'L'v'J',-M'}(\vec{r})\sigma_{v'\mu'}D_{mIM_I}.$$ (59)

In this form, the mean-field operator is expressed through sums of derivative operators standing on both sides of the potentials $\tilde{U}_{n'L'v'J',-M'}$. This form can easily be used in the calculation of the matrix elements, because the left derivative operator can simply be applied onto the left wave-function through the integration by parts. Moreover, potentials $\tilde{U}_{n'L'v'J',-M'}$ are related to the secondary densities by very simple relations (46). However, it turns out that numerical calculations are much faster if the derivatives appear only on one side of the potentials, like it was postulated in Eq. (7), that is,

$$h(\rho) = \sum_{n'L'v'J'}\left[U _{n'L'v'J'}(\vec{r})[D_{n'L'}\sigma_{v'}]_{J'}\right]_0$$

$$= \sum_{n'L'v'J'} \sum_{M'_L\mu'} {\textstyle{\frac{(-1)^{J'-M'}}{\... ...} C^{J'M'}_{L'M'_Lv'\mu'} U_{n'L'v'J',-M'}(\vec{r})\sigma_{v'\mu'}D_{n'L'M'_L}.$$ (60)