text:
# Separable RPA for self-consistent nuclear models

## I Introduction

Self-consistent nuclear models, as e.g. Skyrme-Hartree-Fock (SHF), are widely used for the description of nuclear ground state properties. They do also allow a description of excitation spectra. This is usually done within the random-phase approximation (RPA), which is meanwhile a textbook standard in nuclear physics Row ; Ber94aB , for a few recent applications see Rew ; Col ; Vre . Actually, these dynamical applications are up to now rather limited. The reason is that RPA implies the diagonalization of large matrices whose rank is determined by the size of the particle-hole $`1ph`$ space. The number of $`1ph`$ states grows huge for heavy and/or deformed nuclei. This limits the range of application for fully fledged RPA. The most applications are thus found for spherical nuclei.

The RPA problem becomes much simpler if the residual two-body interaction is reduced to a separable form, i.e.

$`\widehat{V}_{\mathrm{res}}^{}\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k,k^{}=1}{\overset{N_{sep}}{}}}\mu _{kk^{}}\widehat{Z}_k\widehat{Z}_k^{},`$
$`\widehat{Z}_k`$ $`=`$ $`{\displaystyle \underset{ph}{}}p|\widehat{Z}_k|ha_p^{}a_h`$ (1)
where $`\{\widehat{Z}_k\}`$ is a set of hermitian one-body operators. The factorization changes the rank of the RPA matrix from the number of $`1ph`$ configurations to the number of basis operators $`\widehat{Z}_k`$ and this reduces dramatically the computational effort. The success of a separable ansatz depends, of course, on a diligent choice of the operators $`\widehat{Z}_k`$ and the associated strength coefficients $`\mu _{kk^{}}`$.

Factorization of the residual interaction is widely used in nuclear theory but mainly within trivial schemes exploiting one separable term with an intuitive guess for the separable one-body operator $`\widehat{Z}`$. The strength constant $`\mu `$ is usually fitted to reproduce available experimental data (see e.g. Sol ). Obviously, accuracy and predictive power of such simple schemes are limited. Several improvements towards self-consistent schemes have been proposed during the last decades Row -Vor . However, these schemes are not sufficiently general. Some of them are limited to analytic or simple numerical estimates Row -SS , others are not fully self-consistent in the sense that they start from phenomenological single-particle potentials Ne\_PRC -Kubo  or cover only particular effective forces Vor .

In the present paper we propose a general self-consistent separable RPA (SRPA) approach relevant to arbitrary density- and current-dependent functionals. The self-consistent scheme of Ref. LS  is generalized to the case of several separable operators. The operators are chosen to have maxima at different areas of the nucleus. This is crucial for an accurate reproduction of the residual interaction $`V_{\mathrm{res}}`$. A similar scheme has been successfully applied to the Kohn-Sham functional for description of collective oscillations of valence electrons in atomic clusters Ne\_PRA -babst . The nuclear case is more demanding since the SHF functional is much more involved. We will discuss the actual SRPA scheme for the case of SHF and present first successful tests for isoscalar E2 and isovector E1 giant resonances in $`{}_{}{}^{40}\mathrm{Ca}`$ and $`{}_{}{}^{208}\mathrm{Pb}`$ using two typical Skyrme parameterizations.

## II SRPA

### II.1 Separable operators

In connection with nuclear density functional theory, it is preferable to sort the one-body operators $`\widehat{Z}_k`$ according to time parity. We thus decompose the set into hermitian time-even $`\widehat{X}_k`$ and time-odd $`\widehat{Y}_k`$ operators and the corresponding strength matrices $`\kappa _{kk^{}}`$ and $`\eta _{kk^{}}`$. Then the separable Hamiltonian can be written as

$`\widehat{H}_{\mathrm{SRPA}}`$ $`=`$ $`\widehat{h}_0+\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})},`$ (2)
$`\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{kk^{}=1}{\overset{N_{sep}}{}}}(\kappa _{kk^{}}\widehat{X}_k\widehat{X}_k^{}+\eta _{kk^{}}\widehat{Y}_k\widehat{Y}_k^{})`$ (3)
where $`\widehat{h}_0`$ is the mean field Hamiltonian of the ground state, $`\kappa _{kk^{}}`$ and $`\eta _{kk^{}}`$ are symmetric matrices. Note that the operator product is considered to be separable at the Hartree level, i.e. for an arbitrary $`1ph`$ operator $`\widehat{A}`$ we have

$$[\widehat{Z}_k\widehat{Z}_k^{},\widehat{A}]=\widehat{Z}_k[\widehat{Z}_k^{},\widehat{A}]+\widehat{Z}_k^{}[\widehat{Z}_k,\widehat{A}].$$
(4)
Furthermore, it is important to note for the following considerations that only commutators between time-even and time-odd operators have non-vanishing expectation values. That means that $`[\widehat{X}_k,\widehat{X}_k^{}]=[\widehat{Y}_k,\widehat{Y}_k^{}]=0`$ while $`[\widehat{X}_k,\widehat{Y}_k^{}]0`$.

The success of a separable ansatz depends very much on an appropriate choice of the expansion basis $`\widehat{X}_k`$ and $`\widehat{Y}_k`$. Collective modes (vibrational, giant resonances) dominate the spectral strengths and it is intuitively clear that collective operators like the various multipole moments should play the leading role in the expansion. We will denote them as hermitian time-even collective operators $`\widehat{Q}_k`$ (as e.g. a local multipole operator $`\widehat{Q}_{\lambda \mu }(𝐫)r^\lambda (Y_{\lambda \mu }+Y_{\lambda \mu }^{})`$) and time-odd conjugate momenta

$$\widehat{P}_k=i[\widehat{H},\widehat{Q}_k]_{ph}$$
(5)
where the index $`ph`$ means the $`1ph`$ part of the operator. These operators can serve as generators of a collective motion by virtue of the scaling transformation

$$|\mathrm{\Psi }=e^{\widehat{G}}|,\widehat{G}=_k(ip_k\widehat{Q}_kiq_k\widehat{P}_k)$$
(6)
where $`||\mathrm{\Psi }_0`$ is the ground state wave function, $`\widehat{G}`$ is the generator, composed of various coordinate-like $`\widehat{Q}_k`$ and momentum-like $`\widehat{P}_k`$ operators, $`p_k`$ and $`q_k`$ are the corresponding (time-dependent) amplitudes, see e.g. Boh85 ; Bra89 ; ReiBraGen . The dynamical deformation (6) induces changes of the mean-field Hamiltonian $`\widehat{h}`$ which can be expanded in the linear regime for the small $`\widehat{G}`$ as

$$\begin{array}{ccccc}\hfill \widehat{h}& & \widehat{h}_0\hfill & +& \underset{}{[\widehat{h}_0,\widehat{G}]_{ph}+[\widehat{V}_{\mathrm{res}},\widehat{G}]_{ph}}\\ & =& \widehat{h}_0\hfill & +& [\widehat{H},\widehat{G}]_{ph}.\end{array}$$
(7)
The dynamics of small-amplitude vibrations is governed by $`[\widehat{H},\widehat{G}]_{ph}`$ and so we optimize the separable interaction by requiring that

$$[\widehat{H},\widehat{G}]_{ph}=[\widehat{H}^{(\mathrm{sep})},\widehat{G}]_{ph}$$
(8)
as close as possible. Since both Hamiltonians contain $`\widehat{h}_0`$ in the same way, the requirement (8) is essentially

$$[\widehat{V}_{\mathrm{res}},\widehat{G}]_{ph}=[\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})},\widehat{G}]_{ph}.$$
(9)
Applying (9) to generators $`\widehat{Q}_k`$ and $`\widehat{P}_k`$ separately, we have

$`[\widehat{V}_{\mathrm{res}},\widehat{P}_k]_{ph}`$ $`=`$ $`[\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})},\widehat{P}_k]_{ph}`$ (10a)
$`=`$ $`{\displaystyle \underset{\stackrel{~}{k}k^{}}{}}\widehat{X}_{\stackrel{~}{k}}\kappa _{\stackrel{~}{k}k^{}}[\widehat{X}_k^{},\widehat{P}_k],`$
$`[\widehat{V}_{\mathrm{res}},\widehat{Q}_k]_{ph}`$ $`=`$ $`[\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})},\widehat{Q}_k]_{ph}`$ (10b)
$`=`$ $`{\displaystyle \underset{\stackrel{~}{k}k^{}}{}}\widehat{Y}_{\stackrel{~}{k}}\eta _{\stackrel{~}{k}k^{}}[\widehat{Y}_k^{},\widehat{Q}_k].`$
This shows that an optimal choice has to fulfill the condition

$`\widehat{X}_k\left\{[\widehat{V}_{\mathrm{res}},\widehat{P}_k]_{ph}\right\},`$ (11a)
$`\widehat{Y}_k\left\{[\widehat{V}_{\mathrm{res}},\widehat{Q}_k]_{ph}\right\}.`$ (11b)
We are free to assume any linear recombination among the sets $`\{\widehat{X}_k,\widehat{Y}_k\}`$ and thus we may choose

$`\widehat{X}_k`$ $`=`$ $`[\widehat{V}_{\mathrm{res}},\widehat{P}_k]_{ph},`$ (12a)
$`\widehat{Y}_k`$ $`=`$ $`[\widehat{V}_{\mathrm{res}},\widehat{Q}_k]_{ph}.`$ (12b)
Substitution of $`\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})}`$ into (12a) and (12b) has to give the identities. This results in the expressions for the inverse strength constants $`\kappa _{k^{}\overline{k}}^1`$ and $`\eta _{k^{}\overline{k}}^1`$ ($`_{\overline{k}}\kappa _{k\overline{k}}^1\kappa _{\overline{k}k^{}}=\delta _{kk^{}}`$, $`_{\overline{k}}\eta _{k\overline{k}}^1\eta _{\overline{k}k^{}}=\delta _{kk^{}}`$):

$`\kappa _{k^{}\overline{k}}^1`$ $`=`$ $`[\widehat{X}_k^{},\widehat{P}_{\overline{k}}]=[[\widehat{V}_{\mathrm{res}},\widehat{P}_k]_{ph},\widehat{P}_{\overline{k}}],`$ (13a)
$`\eta _{k^{}\overline{k}}^1`$ $`=`$ $`[\widehat{Y}_k^{},\widehat{Q}_{\overline{k}}]=[[\widehat{V}_{\mathrm{res}},\widehat{Q}_k]_{ph},\widehat{Q}_{\overline{k}}].`$ (13b)
The actual calculation of the commutators can be done either in wavefunction representation ReiGam ; rpanucl  or in terms of explicit $`1ph`$ matrix elements. The latter reads

$`[\widehat{X}_k,\widehat{P}_k^{}]=2{\displaystyle \underset{ph}{}}(p|\widehat{X}_k|h)(h|\widehat{P}_k^{}|p),`$ (14a)
$`[\widehat{Y}_k,\widehat{Q}_k^{}]=2{\displaystyle \underset{ph}{}}(p|\widehat{Y}_k|h)(h|\widehat{Q}_k^{}|p).`$ (14b)
Alternatively, the strength constants can be determined from the requirement

$`[\widehat{P}_k,[\widehat{V}_{\mathrm{res}},\widehat{P}_k^{}]]`$ $`=`$ $`[\widehat{P}_k,[\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})},\widehat{P}_k^{}]],`$ (15a)
$`[\widehat{Q}_k,[\widehat{V}_{\mathrm{res}},\widehat{Q}_k^{}]]`$ $`=`$ $`[\widehat{Q}_k,[\widehat{V}_{\mathrm{res}}^{(\mathrm{sep})},\widehat{Q}_k^{}]].`$ (15b)
Substituting (3) and (12a) into (15a) yields

$$[\widehat{P}_k,\widehat{X}_k^{}]=\underset{\stackrel{~}{k}\stackrel{~}{k}^{}}{}[\widehat{P}_k,\widehat{X}_{\stackrel{~}{k}}]\kappa _{\stackrel{~}{k}\stackrel{~}{k}^{}}[\widehat{X}_{\stackrel{~}{k}^{}},\widehat{P}_k^{}]$$
(16)
and finally Eq. (13a). A similar procedure applies to the double commutator (15b). For a physical interpretation of the requirements (15a) and (15b) mind that the similar double commutators but with the full Hamiltonian (instead of the residual interaction) correspond to $`m_3`$ and $`m_1`$ sum rules, respectively, and so represent the spring and inertia parameters ReiBraGen  in the basis of collective generators $`\widehat{Q}_k`$ and $`\widehat{P}_k`$. The condition (13) means that collective moments $`m_3`$ and $`m_1`$ are exactly reproduced.

The basic result of this section are Eqs. (12) and (13) for the separable operators and strength constants, respectively. It is worth noting that these results can be also obtained Maiori  using the concept of nuclear self-consistency between the single-particle potential and density Row . More detailed formulation of the basic double commutators in terms of the SHF functional will be given in section II.5 and appendix B. The links to the vibrating potential model (VPM) Row ; BM  are outlined in appendix A. The choice of the generators $`\widehat{Q}_k`$ and $`\widehat{P}_k`$ is yet open. An intuitive and pragmatic guess will be presented and discussed in sections II.6 and III.

### II.2 Conservation laws

The forms (12) are able to reproduce given conservation laws. Assume a symmetry mode whose generator $`\widehat{P}_{\mathrm{sym}}`$ obeys by definition $`[\widehat{H},\widehat{P}_{\mathrm{sym}}]=0`$. We simply have to include this mode into the set of generators, i.e. $`\widehat{P}_{\mathrm{sym}}\{\widehat{P}_k\}`$ together with its complement $`\widehat{Q}_{\mathrm{sym}}`$ given by $`[\widehat{H},\widehat{Q}_{\mathrm{sym}}]=i\widehat{P}_{\mathrm{sym}}`$. One can prove that the separable Hamiltonian (2) does also fulfill

$$[\widehat{H}_{\mathrm{SRPA}},\widehat{P}_{\mathrm{sym}}]=0.$$
(17)
From (17) we have

$$[\widehat{h}_0,\widehat{P}_{\mathrm{sym}}]=[\widehat{V}_{\mathrm{res}},\widehat{P}_{\mathrm{sym}}]=\widehat{X}_{\mathrm{sym}}.$$

Then

$`[\widehat{H}_{\mathrm{SRPA}},\widehat{P}_{\mathrm{sym}}]`$ $`=`$ $`[\widehat{h}_0,\widehat{P}_{\mathrm{sym}}]+{\displaystyle \underset{k,k^{}}{}}\widehat{X}_k\kappa _{kk^{}}[\widehat{X}_k^{},\widehat{P}_{\mathrm{sym}}]`$
$`=`$ $`\widehat{X}_{\mathrm{sym}}+{\displaystyle \underset{k}{}}\widehat{X}_k\underset{\delta _{k,\mathrm{sym}}}{\underset{}{{\displaystyle \underset{k^{}}{}}\kappa _{kk^{}}\kappa _{k^{}\mathrm{sym}}^1}}`$
$`=`$ $`\widehat{X}_{\mathrm{sym}}+\widehat{X}_{\mathrm{sym}}`$
which obviously yields Eq. (17). Thus all symmetry modes are recovered if properly included in the ansatz.

It is worth noting that in the case of a symmetry mode the operator of the residual interaction is

$$\widehat{X}_{\mathrm{sym}}=[\widehat{h}_0,\widehat{P}_{\mathrm{sym}}],$$
(18)
i.e. can be also presented as the commutator with the single-particle Hamiltonian. In this case, the residual interaction serves to restore the symmetries violated in the mean field. The similar forces were proposed in Pyatov  to restore translational and rotational invariance.

### II.3 Response functions

One way to solve the SRPA equations is to compute directly a desired strength function

$`S_D(\omega )`$ $`=`$ $`{\displaystyle \underset{N}{}}|\mathrm{\Psi }_N|\widehat{D}|\mathrm{\Psi }_0|^2(\delta (\omega \omega _N)+\delta (\omega +\omega _N))`$ (19)
$`=`$ $`{\displaystyle \frac{1}{\pi }}\mathrm{}\left\{[\widehat{D},{\displaystyle \frac{1}{\omega +i\mathrm{\Gamma }}}\widehat{D}]\right\}`$
where $`\widehat{D}`$ is the operator of interest, $`\mathrm{\Psi }_N`$ and $`\omega _N`$ are the wave function and energy of $`N`$-th RPA state, $``$ is the Liouvillian for small amplitude excitations. The latter is a shorthand for

$$\widehat{A}=[\widehat{H},\widehat{A}]_{ph}$$
(20)
where $`\widehat{A}`$ is any $`1ph`$ operator. Similarly, we introduce the unperturbed Liouvillian $`_0`$ and the Liouvillian for the residual interaction $`𝒱_{\mathrm{res}}`$ by

$$_0\widehat{A}=[\widehat{h}_0,\widehat{A}],𝒱_{\mathrm{res}}\widehat{A}=[\widehat{V}_{\mathrm{res}},\widehat{A}]_{ph}.$$

For a compact formulation, we collect the operators $`\widehat{D}`$, $`\{\widehat{X}_k,k=1,\mathrm{},N_{\mathrm{sep}}\}`$ and $`\{\widehat{Y}_k,k=1,\mathrm{},N_{\mathrm{sep}}\}`$ into the set $`\{\widehat{Z}_k,k=0,\mathrm{},2N_{\mathrm{sep}}\}`$ with $`\widehat{X}_0=\widehat{D}`$. The matrices of the strength constants are joined into the matrix $`\mu `$ of the rank $`2N_{\mathrm{sep}}+1`$:

$$\mu =\left(\begin{array}{ccc}1& 0& 0\\ 0& \kappa & 0\\ 0& 0& \eta \end{array}\right).$$
(21)
Then we define generalized response functions

$`R_{kk^{}}(\omega )`$ $`=`$ $`[\widehat{Z}_k,{\displaystyle \frac{1}{\omega +i\mathrm{\Gamma }}}\widehat{Z}_k^{}],`$ (22)
$`R_{kk^{}}^{(0)}(\omega )`$ $`=`$ $`[\widehat{Z}_k,{\displaystyle \frac{1}{\omega +i\mathrm{\Gamma }_0}}\widehat{Z}_k^{}].`$ (23)
Note that the free response function can be written easily in terms of single-particle states $`p`$ and $`h`$ as

$`R_{kk^{}}^{(0)}(\omega )`$ $`=`$ $`{\displaystyle \underset{ph}{}}(p|\widehat{Z}_k|h)(h|\widehat{Z}_k^{}|p)`$ (24)
$`\left\{{\displaystyle \frac{1}{\omega +i\mathrm{\Gamma }\epsilon _{ph}}}{\displaystyle \frac{\mathrm{\Pi }_k\mathrm{\Pi }_k^{}}{\omega +i\mathrm{\Gamma }+\epsilon _{ph}}}\right\}.`$
where $`\epsilon _{ph}`$ is the energy of the 1ph-configuration and $`\mathrm{\Pi }_k`$ is the time-parity of the related operator $`\widehat{Z}_k`$.

A relation between full and free responses is established through the following steps:

$`R_{kk^{}}`$ $`=[\widehat{Z}_k,{\displaystyle \frac{1}{\omega +i\mathrm{\Gamma }_0𝒱_{\mathrm{res}}^{(\mathrm{sep})}}}\widehat{Z}_k^{}]`$
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}[\widehat{Z}_k,{\displaystyle \frac{1}{\omega +i\mathrm{\Gamma }_0}}\left(𝒱_{\mathrm{res}}^{(\mathrm{sep})}{\displaystyle \frac{1}{\omega +i\mathrm{\Gamma }_0}}\right)^n\widehat{Z}_k^{}]`$
$`=`$ $`R_{kk^{}}^{(0)}+{\displaystyle \underset{k_1k_2}{}}R_{kk_1}^{(0)}[\widehat{Z}_{k_1},[𝒱_{\mathrm{res}}^{(\mathrm{sep})},\widehat{Z}_{k_2}]]R_{k_2k^{}}^{(0)}+\mathrm{}`$
$`=`$ $`R_{kk^{}}^{(0)}+{\displaystyle \underset{k_1k_2}{}}R_{kk_1}^{(0)}\mu _{k_1k_2}R_{k_2k^{}}^{(0)}+\mathrm{}.`$
This can be recollected to the final result

$$R_{kk^{}}(\omega )=\left(R^{(0)}(\omega )\frac{1}{1\mu R^{(0)}(\omega )}\right)_{kk^{}}.$$
(26)
Note that this operator relation deals with matrices of the rank $`2(2N_{\mathrm{sep}}+1)`$ (the additional coefficient 2 arises due to isospin). The strength function is finally

$$S_D(\omega )=\frac{1}{\pi }\mathrm{}\left\{R_{00}(\omega )\right\}.$$
(27)
The calculation of the free response function requires some caution with respect to time parity, see Eq. (23). The denominator in $`R_{kk^{}}^{(0)}`$ mixes time parities and one has to disentangle the pieces

$`{\displaystyle \frac{1}{\omega +i\mathrm{\Gamma }_0}}`$ $`=`$ $`_+^{(0)}(\omega )+_{}^{(0)}(\omega ),`$
$`_+^{(0)}(\omega )`$ $`=`$ $`{\displaystyle \frac{_0}{(\omega +i\mathrm{\Gamma })^2_0^2}},`$
$`_{}^{(0)}(\omega )`$ $`=`$ $`{\displaystyle \frac{\omega +i\mathrm{\Gamma }}{(\omega +i\mathrm{\Gamma })^2_0^2}}`$
where $`_+^{(0)}`$ and $`_{}^{(0)}`$ connect operators of the same and opposite parity, respectively. One thus has to use $`_+^{(0)}`$ for $`[\widehat{X},^{(0)}\widehat{X}]`$, $`[\widehat{Y},^{(0)}\widehat{Y}]`$, $`[\widehat{X},^{(0)}\widehat{D}]`$, $`[\widehat{D},^{(0)}\widehat{D}]`$ and $`_{}^{(0)}`$ for $`[\widehat{X},^{(0)}\widehat{Y}]`$, $`[\widehat{D},^{(0)}\widehat{Y}]`$. The mutually other combinations vanish.

### II.4 SRPA as eigenvalue problem

As often done in RPA, one wants to determine the detailed eigenfrequencies and eigenstates of the response problem. The excitation phonon operator $`\widehat{C}_N^{}`$ for the mode $`N`$ is determined from the RPA equations

$$[\widehat{H},\widehat{C}_N^{}]_{ph}=\omega _N\widehat{C}_N^{}.$$
(28)
The standard RPA scheme proceeds to recast that equations into the form of a matrix equation by expanding $`\widehat{C}_N^{}`$ in terms of $`1ph`$ operators

$$\widehat{C}_N^{}=\underset{ph}{}(c_{N,ph}^{(+)}\widehat{a}_p^{}\widehat{a}_h^{}c_{N,ph}^{()}\widehat{a}_h^{}\widehat{a}_p^{}).$$
(29)
The involved matrices can grow huge in view of the large $`1ph`$ spaces required for sufficient convergence of the result. The eigenvalue problem is dramatically simplified by inserting the separable Hamiltonian (2). For compact notation, we recombine again the sets $`\{\widehat{X}_k,\widehat{Y}_k\}\{\widehat{Z}_k,k=1,\mathrm{},2N_{\mathrm{sep}}\}`$ and introduce the corresponding super-matrix of strengths $`(\kappa _{kk^{}},\eta _{kk^{}})\mu _{kk^{}}`$. The RPA equations (28) then become

$`c_{N,ph}^{(+)}`$ $`=`$ $`{\displaystyle \underset{kk^{}}{}}\mu _{kk^{}}{\displaystyle \frac{(p|\widehat{Z}_k|h)𝒞_{N,k^{}}}{\omega _N\epsilon _{ph}}},`$ (30a)
$`c_{N,ph}^{()}`$ $`=`$ $`{\displaystyle \underset{kk^{}}{}}\mu _{kk^{}}{\displaystyle \frac{(h|\widehat{Z}_k|p)𝒞_{N,k^{}}}{\omega _N+\epsilon _{ph}}}`$ (30b)
with

$`𝒞_{N,k}`$ $`=`$ $`[\widehat{Z}_k,\widehat{C}_N^{}]`$
$`=`$ $`{\displaystyle \underset{ph}{}}\left\{(h|\widehat{Z}_k|p)c_{N,ph}^{(+)}+(p|\widehat{Z}_k|h)c_{N,ph}^{()}\right\}.`$
Eqs. (30a)-(30b) reduce the impressive number of unknowns $`c_{N,ph}^{(\pm )}`$ to the much smaller set of unknonws $`𝒞_{N,k}`$. Such reduction (which is possible due to factorization of the residual interaction) yields an effective RPA matrix with much smaller rank.

As a next step, it is then convenient to transform $`𝒞_{N,k}`$ to unknowns

$$\overline{𝒞}_{N,k}=2\underset{k^{}}{}\mu _{kk^{}}𝒞_{N,k^{}}$$
(32)
which are directly connected with the amplitudes $`p_k(t)`$ and $`q_k(t)`$ of the scaling transformation (6). For harmonic oscillations

$$q_k(t)=\overline{q}_kcos(\omega t),p_k(t)=\overline{p}_ksin(\omega t)$$
(33)
and $`\overline{𝒞}_{N,k}=\overline{q}_k`$ and $`\overline{p}_k`$ for time-even and time-odd cases, respectively. The unknowns (32) regulate self-consistently contributions of the scaling operators $`\widehat{Q}_k`$ and $`\widehat{P}_k`$ into N-th RPA state. Every RPA state has it own set $`\overline{𝒞}_{N,k}`$ ($`\overline{q}_k`$ and $`\overline{p}_k`$) amplitudes.

Finally, after standard algebraic steps, we obtain the SRPA equation for $`\overline{𝒞}_{N,k}`$

$$\underset{k^{}}{}d_{kk^{}}(\omega _N)\overline{𝒞}_{N,k^{}}=0$$
(34)
with

$$d_{kk^{}}(\omega _N)=2\underset{ph}{}\frac{(p|\widehat{Z}_k|h)(h|\widehat{Z}_k^{}|p)\epsilon _{ph}}{\epsilon _{ph}^2\omega _N^2}\mu _{kk^{}}^1.$$
(35)
Eq. (34) has non-trivial solutions if

$$\text{det}\left\{d_{kk^{}}(\omega _N)\right\}=0$$
(36)
which yields the RPA spectrum $`\omega _N`$. Eq. (34) is basically a matrix of the rank $`4N_{\mathrm{sep}}`$ (after duplication of the rank because of the isospin) which is much smaller than the necessary number of $`1ph`$ states. Typically, we can deal with $`N_{\mathrm{sep}}=14`$ (see section III) as opposed to the rank of a conventional RPA matrix which often amounts to about $`10^310^6`$.

The normalization condition for RPA states is

$$(\widehat{C}_N^{},\widehat{C}_N^{})=1$$
(37)
which amounts to

$`{\displaystyle \underset{ph}{}}(|c_{N,ph}^{(+)}|^2|c_{N,ph}^{()}|^2)`$
$`=`$ $`{\displaystyle \underset{kk^{}}{}}\overline{𝒞}_{N,k}\overline{𝒞}_{N,k^{}}{\displaystyle \underset{ph}{}}{\displaystyle \frac{<p|\widehat{Z}_k|h><p|\widehat{Z}_k^{}|h>\epsilon _{ph}\omega _N}{(\epsilon _{ph}^2\omega _j^2)^2}}`$
$`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{kk^{}}{}}\overline{𝒞}_{N,k}\overline{𝒞}_{N,k^{}}{\displaystyle \frac{}{\omega _N}}d_{kk^{}}(\omega _N)=2.`$
The last line of (II.4) represents the normalization condition in terms of derivatives of the RPA matrix elements.

The reduced probability of $`E\lambda `$-transition from the ground state to the RPA state $`\omega _N`$ can be written as

$$B(E\lambda ,gr\omega _N)=\frac{_{kk^{}}D_{kk^{}}(\omega _N)A_k(\omega _N)A_k^{}(\omega _N)}{\frac{}{\omega _N}\text{det}\left\{d_{kk^{}}(\omega _N)\right\}}$$
(39)
where $`D_{kk^{}}(\omega _N)`$ is the algebraic supplement of the matrix element $`d_{kk^{}}(\omega _N)`$,

$$A_k(\omega _N)=\underset{ph}{}\frac{\epsilon _{ph}(p|\widehat{Z}_k|h)(h|\widehat{f}_\lambda |p)}{\epsilon _{ph}^2\omega _N^2}$$
(40)
and $`(h|\widehat{f}_\lambda |p)`$ is the 1ph matrix element of $`E\lambda `$-transition.

### II.5 Connection with density functional theory

Self-consistent models are often formulated in terms of a local energy-density functional rather than through an explicit Hamiltonian (see, e.g., the nuclear Skyrme-Hartree-Fock (SHF) method SHFrev  or the Kohn-Sham scheme in electron systems KohSha ; DreGro ). So we will work out the above SRPA scheme in terms of a given energy functional

$$E=E(\{J_\alpha (𝐫)\})$$

where

$$J_\alpha (𝐫)=\mathrm{\Psi }|\widehat{J}_\alpha (𝐫)|\mathrm{\Psi }$$

constitute a set of local densities and currents associated with the functional. A simple example is the situation in atoms and atomic clusters where we deal only with the local density $`\rho (𝐫)=_n|\phi _n(𝐫)|^2`$ of valence electrons. The nuclear case is much more involved. The typical set for Skyrme functional includes, for both protons and neutrons, the local density $`\rho `$, kinetic density $`\tau `$, current $`𝐣`$, spin-orbit density $`𝐉`$, and spin density $`\sigma `$ i.e.

$$J_\alpha \{\rho _p,\rho _n,\tau _p,\tau _n,𝐉_p,𝐉_n,𝐣_p,𝐣_n,\sigma _p,\sigma _n\},$$
(41)
see in Appendix B the detailed expressions.

The mean field Hamiltonian corresponding to the given energy functional is then determined as

$$\widehat{h}(𝐫)=𝑑𝐫\underset{\alpha }{}\widehat{J}_\alpha (𝐫)\frac{\delta E}{\delta J_\alpha (𝐫)}.$$
(42)
In the linear regime of small-amplitude deformation generated by a $`1ph`$ operator $`\widehat{G}`$, see Eq. (6), the densities can be decomposed as

$$J_\alpha (𝐫,t)=\overline{J}_\alpha (𝐫)+J_\alpha ^G(𝐫,t)$$
(43)
where $`\overline{J}_\alpha (𝐫)`$ is the static ground-state density and

$`J_\alpha ^G(𝐫,t)`$ $`=`$ $`\mathrm{\Psi }(t)|\widehat{J}_\alpha (𝐫)|\mathrm{\Psi }(t)\mathrm{\Psi }_0|\widehat{J}_\alpha (𝐫)|\mathrm{\Psi }_0`$ (44)
$`=`$ $`[\widehat{J}_\alpha (𝐫),\widehat{G}(𝐫,t)]`$
accounts for the small change through the deformation. Inserting (43) into the mean-field Hamiltonian (42) and expanding the latter in orders of $`\widehat{G}`$, we get, in the first order, the response Hamiltonian $`[\widehat{V}_{\mathrm{res}},\widehat{G}]_{ph}^{}`$ in terms of the density functional

$$[\widehat{V}_{\mathrm{res}},\widehat{G}]_{ph}^{}=\underset{\alpha \beta }{}𝑑𝐫𝑑𝐫^{}\widehat{J}_\alpha (𝐫)\frac{\delta ^2E}{\delta J_\alpha (𝐫)\delta J_\beta (𝐫^{})}J_\beta ^G(𝐫^{},t).$$
(45)
We now decompose the scaling state as in Eq. (6), which yields the analogous decomposition for the variation of density

$$J_\alpha ^G(𝐫,t)=i\underset{k}{}\left(p_k(t)J_\alpha ^{Q_k}(𝐫)q_k(t)J_\alpha ^{P_k}(𝐫)\right)$$
(46)
where, similar to Eq. (44),

$`J_\alpha ^{Q_k}(𝐫)`$ $`=`$ $`[\widehat{J}_\alpha (𝐫),\widehat{Q}_k(𝐫)],`$ (47a)
$`J_\alpha ^{P_k}(𝐫)`$ $`=`$ $`[\widehat{J}_\alpha (𝐫),\widehat{P}_k(𝐫)].`$ (47b)
Every step proceeds much similar as for the general $`\widehat{G}`$ because we are in the linear regime. Finally, we can specify the separable operators (12a) and (12b) as

$`\widehat{X}_k(𝐫)`$ $`=`$ $`{\displaystyle \underset{\alpha \beta }{}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\widehat{J}_\alpha (𝐫)\frac{\delta ^2E}{\delta J_\alpha (𝐫)\delta J_\beta (𝐫^{})}J_\beta ^{P_k}(𝐫^{})},`$ (48a)
$`\widehat{Y}_k(𝐫)`$ $`=`$ $`{\displaystyle \underset{\alpha \beta }{}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\widehat{J}_\alpha (𝐫)\frac{\delta ^2E}{\delta J_\alpha (𝐫)\delta J_\beta (𝐫^{})}J_\beta ^{Q_k}(𝐫^{})}.`$ (48b)
and the strength matrices (13a) and (13b) as

$`\kappa _{kk^{}}^1`$ $`=`$ $`{\displaystyle \underset{\alpha \beta }{}}{\displaystyle 𝑑𝐫𝑑𝐫^{}J_\alpha ^{P_k}(𝐫)\frac{\delta ^2E}{\delta J_\alpha (𝐫)\delta J_\beta (𝐫^{})}J_\beta ^{P_k^{}}(𝐫^{})},`$ (49a)
$`\eta _{kk^{}}^1`$ $`=`$ $`{\displaystyle \underset{\alpha \beta }{}}{\displaystyle 𝑑𝐫𝑑𝐫^{}J_\alpha ^{Q_k}(𝐫)\frac{\delta ^2E}{\delta J_\alpha (𝐫)\delta J_\beta (𝐫^{})}J_\beta ^{Q_k^{}}(𝐫^{})}.`$ (49b)
The latter equations demonstrate the symmetry of the strength matrices $`\kappa _{kk^{}}`$ and $`\eta _{kk^{}}`$. Note the subtle difference between separable operators and strengths. The double integrals look much the same. However, in Eqs. (48a) and (48b), they contain the current operators $`\widehat{J}_\alpha (𝐫)`$ and yield one-body operators as a result of the integration. Instead, Eqs. (49a) and (49b) do not include operators and so yield strength coefficients as c-numbers. For more details see Appendix B.

### II.6 Choice of collective generators

As was mentioned above, the proper choice of generating operators $`\widehat{Q}_k(𝐫)`$ is crucial to achieve good convergence of the separable expansion (3) with a minimal number of separable operators. The choice is inspired by both physical and computational arguments. It should be simple and universal in the sense that it can be applied equally well to all modes and excitation channels. The main idea is that the generating operators should explore different spatial regions of the nucleus, the surface as well as the interior. This suggests that the leading scaling operator should have the form of the applied external field in the long-wave approximation, for example, $`\widehat{Q}_{k=1}^{\lambda \mu }(𝐫)=r^\lambda (Y_{\lambda \mu }(\mathrm{\Omega })+\text{h.c.)}`$. Such a choice results in separable operators (48a) and (48b) which are most sensitive to the surface of the system. The first excitation modes of nuclear collective motion tend to vibrate predominantly in the surface region. As a result, already with this one separable term we obtain a quite good description of giant resonances. The detailed distributions depends on a subtle interplay of surface and volume vibrations. This can be resolved by taking into account an interior of the nucleus. For this aim, the radial parts with larger powers $`r^{\lambda +p}Y_{\lambda \mu }`$ and spherical Bessel functions can be used, much similar as in the local RPA ReiBraGen ; rpanucl . This results in the shift of the maxima of the operators (48a) and (48b) to the interior. Exploring different conceivable combinations, we have found a most efficient set of the generators as

$$\widehat{Q}_k(𝐫)=R_k(r)(Y_{\lambda \mu }(\mathrm{\Omega })+\text{h.c.}),\widehat{P}_k=i[\widehat{H},\widehat{Q}_k]$$
(50)
with

$$R_k(r)=\{\begin{array}{cc}r^\lambda ,\hfill & k=1\hfill \\ j_\lambda (q_\lambda ^kr),\hfill & k=2,3,4\hfill \end{array}$$
(51)
$$q_\lambda ^k=a_k\frac{z_\lambda }{R_{\mathrm{diff}}},a_2=0.6,a_3=0.9,a_4=1.2$$

where $`R_{\mathrm{diff}}`$ is the diffraction radius of the actual nucleus and $`z_\lambda `$ is the first root in $`j_\lambda (z_\lambda )=0`$. The separable term with $`k=1`$ is mainly localized at the nuclear surface while the next three terms are localized more and more in the interior. This simple set seems to be a best compromise for the description of nuclear giant resonances in light and heavy nuclei.

One may argue that one needs to explore not only different spatial regions but also the different types of the operators, associated with the variety of densities and currents (e.g. kinetic, spin-orbit). We have checked this. And indeed, this allows a perfect fine tuning of the spectra. On the other hand, it adds more separable terms. The set (50)-(51) with few and purely local generators is in our opinion the best compromise between quality of the results and the expense of the calculations.

## III Results and discussion

We test SRPA for Skyrme forces by comparison with standard full RPA calculations. As the test cases we consider the isovector dipole and isoscalar quadrupole resonances in the two doubly magic nuclei $`{}_{}{}^{40}\mathrm{Ca}`$ and $`{}_{}{}^{208}\mathrm{Pb}`$. Most tests are performed with widely used Skyrme force SkM\* SkM . The convergence of the SRPA expansion is found to be only slightly dependent on the actual force, which is demonstrated for the case of the recent parameterization SkI3 FloRei . The full RPA calculations are done in a large basis which covers the $`1ph`$ excitations from all hole states to all particle states up to 30 MeV above the Fermi energy. SRPA easily allows to extend this basis but we restrict the $`1ph`$ summations in SRPA to precisely the same active space as used in full RPA to have one-to-one comparison. The comparison is done in terms of the strength function $`𝒮_D(\omega )`$ using a smoothing parameter of $`\mathrm{\Gamma }=1\mathrm{MeV}`$. Specifically, the Lorentz weight function was used. The value of the smoothing is quite realistic because experimental spectra are typically washed out by similar widths.

Figure 1 compares SRPA and full RPA results for SkM\* forces. Three stages of SRPA expansion with increasing number of separable terms are shown: only one term $`k=1`$, two terms $`k=1,4`$, and all four terms according to the labeling (51). Mind that we use in the second set the operator with $`a_4=1.2`$ since then the forces cover the extremes of both surface (k=1) and volume (k=4) sensitivity. The first stage with only $`k=1`$ provides already very good description of the isoscalar E2 resonance in both nuclei. Even the low-energy quadrupole mode in $`{}_{}{}^{208}\mathrm{Pb}`$ is well reproduced. The influence of further terms is very small. The isovector E1 resonance seems to be more demanding. The first stage correctly describes the average position and width of the resonance but cannot fully cope with the detailed fragmentation pattern. Here one sees clearly the systematic improvement achieved by the further terms. The set with four terms provides a satisfying agreement. Further improvement is possible by adding more terms but the extra effort seems to be unnecessary. The agreement achieved is already more than sufficient for any practical calculations. One may even be content with a smaller sets for large scale exploratory studies.

Although the smoothing with 1 MeV provides the typical realistic shapes of strength distributions, it is worth having a quick glance at spectra with higher resolution. This is done in Fig. 2 drawn with smoothing $`\mathrm{\Gamma }=`$0.1 MeV. The figure provides the detailed comparison of SRPA and exact RPA results for the most complicated cases of the isovector E1 resonance. Though we have not a perfect coincidence (which may be partly excused by numerical noise in the full RPA calculations), the results are indeed extremely close. This confirms the reliability of SRPA even at this finer scale.

Figure 3 demonstrates the quality of SRPA for another Skyrme force, namely the parameterization SkI3. It shows the toughest test case, the isovector E1 resonance in $`{}_{}{}^{208}\mathrm{Pb}`$. The full set k=1-4 is used. Again we see a satisfying reproduction of this ”worst case”. The other test cases (not shown here) perform even better. We also checked a few other forces and always found the same quality.

In fact, the spectrum from SkI3 looks much more realistic than that obtained from SkM\* (see lower right panel in Fig. 1. The pronounced right shoulder in SkM\* is, to a large extend, an artefact from the restricted $`1ph`$ space used in these calculations. SRPA offers the opportunity to remove that restriction and to perform calculations in the full space of excitations, limited only by numerical resolution of the underlying coordinate space grid (with particle energies at least up to 600 MeV). Figure 4 compares SRPA results in the restricted and unrestricted spaces. It is obvious that the consistent inclusion of $`1ph`$ states above 30 MeV serves to diminish the unnaturally high shoulder at 16 MeV. This is a clear hint that RPA calculations converge slowly with the size of phase space. The SRPA allows to go much further in that respect.

Altogether, the tests show that the proposed set of four operators suffices for a satisfying reproduction of full RPA. The number ”four” has to be taken with a grain of salt. We have four operators of the type $`\widehat{Q}_k`$ for protons and another four for neutrons. The set is complemented by the conjugate operators $`\widehat{P}_k=i[\widehat{H},\widehat{Q}_k]`$. This amounts to 16 operators in total and the computation of the strength function requires inversion of $`17\times 17`$ matrices. The full RPA calculations, on the other hand, used matrices up to a rank of 400 and more.

## IV Conclusions

A novel self-consistent separable approach to the random phase approximation (SRPA) is proposed. SRPA employs a systematic expansion of the exact residual interaction into a sum of weighted separable terms (i.e. products of one-body operators evaluated at a Hartree level). For both, the weights and the separable operators, compact analytical expressions in terms of collective generators are derived in a self-consistent manner. The form of the separable residual interaction is optimized by the condition that the vibrating mean field as generated from collective flow reproduces the exact result. SRPA is formulated in a general way such that it can be applied to arbitrary energy functionals depending on various local densities and currents. As open point remains the choice of the relevant collective generators. Here we exploit experience gained in nuclear fluid dynamical models (local RPA) and use an efficient mix of surface and volume modes.

The general SRPA scheme is specified for the particular case of the Skyrme energy functional. The convergence of the separable expansion and performance of SRPA for giant resonances in light and heavy nuclei is investigated. The results are very encouraging. SRPA proves to be a very good approximation to exact RPA and it it requires much less numerical effort. In particular, it allows to work with huge particle-hole spaces, giving thus the chance for unrestricted RPA calculations, even in nuclei with broken symmetries

## Acknowledgments

This work was partly supported by the Czech grant agency under the contract No. 202/02/0939 (J.K.). V.O.N. thanks the support of RFBR (00-02-17194), Heisenberg-Landau (Germany-BLTP JINR), DFG (436RUS17/102/01) and Votruba-Blokhintcev (Czech Republic -BLTP JINR) grants.

## Appendix A Comparison with VPM

It is interesting to compare SRPA with the Vibrating Potential Model (VPM) Row ; BM  whose different versions are widely used in nuclear structure theory (see, for example, LS ; SS ; Ne\_PRC ; Kubo ; Aberg ). VPM deals with the mean field

$$\widehat{h}_0=T+U_0$$
(52)
where $`U_0`$ is purely local. Moreover, $`U_0`$ is often a density-independent phenomenological potential (oscillator, Nilsson, Woods-Saxon, …). The vibrations of the system with the mean field (52) are characterized by the collective deformations $`q_k`$ and in the framework of VPM the vibrating potential is written as

$$U=U_0+\delta U,\delta U=\underset{k}{}\frac{U}{q_k}\delta q_k=i\underset{k}{}X_k\delta q_k$$
(53)
where $`X_k=i\frac{U}{q_k}`$ are single-particle operators in the r-representation, participating in VPM separable residual interactions (see e.g. Aberg ). Usually the collective momenta $`P_k`$ are introduced

$$P_k=i\frac{}{q_k}.$$
(54)
This allows to write

$`X_k=i{\displaystyle \frac{U}{q_k}}=[U,P_k]`$
$`\widehat{X}_k={\displaystyle 𝑑𝐫\widehat{\rho }(𝐫)[\widehat{U},\widehat{P}_k]}{\displaystyle 𝑑𝐫\widehat{\rho }(𝐫)<[\widehat{U},\widehat{P}_k]>}`$ (55)
where $`\widehat{X}_k`$ is the operator corresponding to $`X_k`$ and acting in the Fock space. VPM operators (A) should be compared with the SRPA operators (12a) or (48a). For the first sight Eq. (A) contradicts to (12a). However, there are, in fact, more similarities than differences. The equivalence of (A) and (48a) can be demonstrated in the case of the simple energy functional

$$E=\frac{1}{2}t_0𝑑𝐫\rho (𝐫)^2.$$
(56)
This functional gives the following mean field potential

$$U=t_0\rho .$$
(57)
VPM expression (A) gives

$$\widehat{X}_k=t_0𝑑𝐫\widehat{\rho }(𝐫)<[\widehat{\rho }(𝐫),\widehat{P}_k]>,$$
(58)
while SRPA expression (48a) provides

$`\widehat{X}_k`$ $`=`$ $`{\displaystyle 𝑑𝐫𝑑𝐫^{}\widehat{\rho }(𝐫)\frac{\delta ^2E}{\delta \rho (𝐫)\delta \rho (𝐫^{})}<[\widehat{\rho }(𝐫^{}),\widehat{P}_k]>}`$ (59)
$`=`$ $`t_0{\displaystyle 𝑑𝐫\widehat{\rho }(𝐫)<[\widehat{\rho }(𝐫),\widehat{P}_k]>}.`$
Thus there is no difference between VPM and SRPA for simple two-body forces following from the functional (56). Differences come up with non-linear terms in the mean-field potential. From this point of view SRPA seems to be more general.

## Appendix B Details of SHF

A most widely used energy functional for nuclear structure calculations is the Skyrme functional. It was originally proposed in skyrme  and made its breakthrough with the first fine-tuning of VauBri . Since then it has been steadily further developed and constitutes today a very reliable model for nuclear structure and excitations. The functional involves the set of local densities and currents (41) which read in general

$`\rho _q(𝐫)`$ $`=`$ $`{\displaystyle \underset{nq}{}}\phi _n^{}(𝐫)\phi _n^{}(𝐫),`$
$`\tau _q(𝐫)`$ $`=`$ $`{\displaystyle \underset{nq}{}}\phi _n^{}(𝐫)\phi _n^{}(𝐫),`$
$`𝐉_q(𝐫)`$ $`=`$ $`i{\displaystyle \underset{nq}{}}\phi _n^{}(𝐫)(\times \widehat{\sigma })\phi _n^{}(𝐫),`$
$`𝐣_q(𝐫)`$ $`=`$ $`{\displaystyle \frac{i}{2}}{\displaystyle \underset{nq}{}}\left[\phi _n^{}(𝐫)\phi _n^{}(𝐫)\phi _n^{}(𝐫)\phi _n^{}(𝐫)\right],`$
$`\sigma _q(𝐫)`$ $`=`$ $`{\displaystyle \underset{nq}{}}\phi _n^{}(𝐫)\widehat{\sigma }\phi _n^{}(𝐫).`$
The associated operators are

$`\widehat{\rho }_q(𝐫)`$ $`=`$ $`\widehat{\mathrm{\Pi }}_q\delta (\widehat{𝐫}^{}𝐫),`$
$`\widehat{\tau }_q(𝐫)`$ $`=`$ $`\widehat{\mathrm{\Pi }}_q\widehat{\stackrel{}{𝐩}}\delta (\widehat{𝐫}^{}𝐫)\widehat{𝐩},`$
$`\widehat{𝐉}_q(𝐫)`$ $`=`$ $`\widehat{\mathrm{\Pi }}_q\delta (\widehat{𝐫}^{}𝐫)\widehat{𝐩}\times \widehat{\sigma },`$
$`\widehat{𝐣}_q(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\widehat{\mathrm{\Pi }}_q\{\widehat{𝐩},\delta (\widehat{𝐫}^{}𝐫)\},`$
$`\widehat{\sigma }_q(𝐫)`$ $`=`$ $`\widehat{\mathrm{\Pi }}_q\delta (\widehat{𝐫}^{}𝐫)\widehat{\sigma }`$
where $`\widehat{\mathrm{\Pi }}_q`$ is the isospin projector to $`q\{p,n\}`$.

We use here the Skyrme functional in the form rpanucl

$`E`$ $`=`$ $`{\displaystyle 𝑑𝐫\left(_{\mathrm{kin}}+_{\mathrm{Sk}}(\rho _q,\tau _q,\sigma _q,𝐣_q,𝐉_q)\right)}`$ (60)
$`+E_\mathrm{C}(\rho _p)E_{\mathrm{cm}},`$
$`_{\mathrm{kin}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\tau ,`$ (61)
$`E_\mathrm{C}`$ $`=`$ $`{\displaystyle \frac{e^2}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\rho _p(𝐫)\frac{1}{|𝐫𝐫^{}|}\rho _p(𝐫^{})}`$ (62)
$`{\displaystyle \frac{3}{4}}e^2({\displaystyle \frac{3}{\pi }})^{\frac{1}{3}}{\displaystyle d^3r[\rho _p(𝐫)]^{\frac{4}{3}}},`$
$`_{\mathrm{Sk}}`$ $`=`$ $`+{\displaystyle \frac{b_0}{2}}\rho ^2{\displaystyle \frac{b_0^{}}{2}}{\displaystyle \underset{q}{}}\rho _q^2+{\displaystyle \frac{b_3}{3}}\rho ^{\alpha +2}{\displaystyle \frac{b_3^{}}{3}}\rho ^\alpha {\displaystyle \underset{q}{}}\rho _q^2+b_1(\rho \tau 𝐣^2)b_1^{}{\displaystyle \underset{q}{}}(\rho _q\tau _q𝐣_q^2)`$ (63)
$`{\displaystyle \frac{b_2}{2}}\rho \mathrm{\Delta }\rho +{\displaystyle \frac{b_2^{}}{2}}{\displaystyle \underset{q}{}}\rho _q\mathrm{\Delta }\rho _qb_4\left(\rho 𝐉+\sigma (\times 𝐣)\right)`$
$`b_4^{}{\displaystyle \underset{q}{}}[\rho _q(𝐉_q)+\sigma _q(\times 𝐣_q)]\stackrel{~}{b}_4𝐉^2\stackrel{~}{b}_4^{}{\displaystyle \underset{q}{}}𝐉_q^2.`$
The sum $`_q`$ runs over proton or neutron species. A density without index means total density, e.g. $`\rho =\rho _p+\rho _n`$ and similarly for the others. The parameters $`b_i`$ and $`b_i^{}`$ used in the above definition are chosen to give a most compact formulation of the energy functional, the corresponding mean-field Hamiltonian and residual interaction . They are related to the standard Skyrme parameters $`t_i`$ and $`x_i`$ VauBri ; SHFrev  by

$$\begin{array}{ccc}\hfill b_0& =& t_0(1+\frac{1}{2}x_0),\hfill \\ \hfill b_0^{}& =& t_0(\frac{1}{2}+x_0),\hfill \\ \hfill b_1& =& \frac{1}{4}\left[t_1(1+\frac{1}{2}x_1)+t_2(1+\frac{1}{2}x_2)\right],\hfill \\ \hfill b_1^{}& =& \frac{1}{4}\left[t_1(\frac{1}{2}+x_1)t_2(\frac{1}{2}+x_2)\right],\hfill \\ \hfill b_2& =& \frac{1}{8}\left[3t_1(1+\frac{1}{2}x_1)t_2(1+\frac{1}{2}x_2)\right],\hfill \\ \hfill b_2^{}& =& \frac{1}{8}\left[3t_1(\frac{1}{2}+x_1)+t_2(\frac{1}{2}+x_2)\right],\hfill \\ \hfill b_3& =& \frac{1}{4}t_3(1+\frac{1}{2}x_3),\hfill \\ \hfill b_3^{}& =& \frac{1}{4}t_3(\frac{1}{2}+x_3),\hfill \\ \hfill b_4& =& b_4^{}=\frac{1}{2}t_4.\hfill \end{array}$$
(64)
Various versions exist for the spin-orbit term. The $`𝐉^2`$ term is sometimes considered, sometimes not. Thus one has the choice

$$\begin{array}{ccc}\hfill \stackrel{~}{b}_4& =& \frac{1}{16}(t_1x_1+t_2x_2),\hfill \\ \hfill \stackrel{~}{b}_4^{}& =& +\frac{1}{16}(t_1t_2)\hfill \end{array}$$
(65)
or

$$\stackrel{~}{b}_4=0,\stackrel{~}{b}_4^{}=0.$$
(66)
The $`\rho 𝐉`$ term is always taken into account. Two options exist here. Conventional Skyrme forces use $`b_4^{}=b_4`$ while a recent extension allows for more flexible isospin dependence considering $`b_4`$ and $`b_4^{}`$ as independent parameters FloRei . The center-of-mass correction is $`E_{\mathrm{cm}}=\widehat{P}_{\mathrm{cm}}^2/2mA`$, often approximated by its diagonal terms $`E_{\mathrm{cm}}=_i\widehat{p}_i^2/2mA`$. The full correction is usually subtracted a posteriori and the approximate correction is expectation value of a one-body operator. In any case, there is no contribution to the residual interaction for RPA.

First variation with respect to $`J_\alpha `$ yields the mean-field Hamiltonian. It can be found in many publications, e.g. rpanucl . The step of interest is the second variation yielding the residual interaction. We obtain here

$`{\displaystyle \frac{\delta ^2E}{\delta \rho _q(𝐫)\delta \rho _q^{}(𝐫^{})}}`$ $`=`$ $`\left(b_0+b_3{\displaystyle \frac{(\alpha +1)(\alpha +2)}{3}}\rho ^\alpha b_3^{}{\displaystyle \frac{\alpha (\alpha 1)}{3}}\rho ^{\alpha 2}{\displaystyle \underset{q}{}}\rho _q^2b_3^{}{\displaystyle \frac{2\alpha }{3}}\rho ^{\alpha 1}b_2\mathrm{\Delta }\right)\delta (𝐫𝐫^{})`$
$`\left(b_0^{}+b_3^{}{\displaystyle \frac{2}{3}}\rho ^\alpha +b_3^{}{\displaystyle \frac{2\alpha }{3}}\rho ^{\alpha 1}\rho _q+b_2^{}\mathrm{\Delta }\right)\delta _{qq^{}}\delta (𝐫𝐫^{})`$
$`+\delta _{qp}\delta _{pq^{}}\left({\displaystyle \frac{e^2}{|𝐫𝐫^{}|}}{\displaystyle \frac{1}{3}}e^2({\displaystyle \frac{3}{\pi }})^{\frac{1}{3}}\rho _p^{2/3}\delta (𝐫𝐫^{})\right),`$
$`{\displaystyle \frac{\delta ^2E}{\delta \rho _q(𝐫)\delta \tau _q^{}(𝐫^{})}}`$ $`=`$ $`b_1\delta (𝐫𝐫^{})b_1^{}\delta (𝐫𝐫^{})\delta _{qq^{}},`$
$`{\displaystyle \frac{\delta ^2E}{\delta \rho _q(𝐫)\delta 𝐉_q^{}(𝐫^{})}}`$ $`=`$ $`\left(b_4b_4^{}\delta _{qq^{}}\right)\delta (𝐫𝐫^{}),`$
$`{\displaystyle \frac{\delta ^2E}{\delta 𝐉_{\nu q}(𝐫)\delta 𝐉_{\nu ^{}q^{}}(𝐫^{})}}`$ $`=`$ $`\left(\stackrel{~}{b}_4+\stackrel{~}{b}_4^{}\delta _{qq^{}}\right)\delta (𝐫𝐫^{})\delta _{\nu \nu ^{}},`$
$`{\displaystyle \frac{\delta ^2E}{\delta j_{\nu q}(𝐫)\delta j_{\nu ^{}q^{}}(𝐫^{})}}`$ $`=`$ $`\left(2b_12b_1^{}\delta _{qq^{}}\right)\delta (𝐫𝐫^{})\delta _{\nu \nu ^{}},`$
$`{\displaystyle \frac{\delta ^2E}{\delta j_{\nu q}(𝐫)\delta \sigma _{\nu ^{}q^{}}(𝐫^{})}}`$ $`=`$ $`\left(b_4+b_4^{}\delta _{qq^{}}\right)_\mu \times \delta (𝐫𝐫^{})ϵ_{\nu \nu ^{}\mu }`$
where $`\nu ,\nu ^{}\{x,y,z\}`$ labels the spatial vector components and $`ϵ_{\nu \nu ^{}\mu }`$ is the antisymmetric tensor. All combinations of derivatives not listed above vanish.

It is instructive to write the response currents in detail. The non-vanishing contributions are

$`\rho _q^{P_k}(𝐫)`$ $`=`$ $`<[\widehat{\rho }_q,\widehat{P}_k]>`$
$`=`$ $`{\displaystyle \underset{nq}{}}\left(\phi _n^{}\widehat{P}_k\phi _n^{}\widehat{P}_k\phi _n^{}\phi _n^{}\right),`$
$`\tau _q^{P_k}(𝐫)`$ $`=`$ $`<[\widehat{\tau }_q,\widehat{P}_k]>`$
$`=`$ $`{\displaystyle \underset{nq}{}}\left(\phi _n^{}\widehat{P}_k\phi _n^{}\widehat{P}_k\phi _n^{}\phi _n^{}\right),`$
$`𝐉_q^{P_k}(𝐫)`$ $`=`$ $`<[\widehat{𝐉}_q,\widehat{P}_k]>`$
$`=`$ $`i{\displaystyle \underset{nq}{}}\left(\phi _n^{}(\times \widehat{\sigma })\widehat{P}_k\phi _n^{}\widehat{P}_k\phi _n^{}(\times \widehat{\sigma })\phi _n^{}\right),`$
$`𝐣_q^{Q_k}(𝐫)`$ $`=`$ $`<[\widehat{𝐣}_q,\widehat{Q}_k]>`$
$`=`$ $`\mathrm{}\left\{{\displaystyle \underset{nq}{}}\left(\phi _n^{}\widehat{Q}_k\phi _n^{}\widehat{Q}_k\phi _n^{}\phi _n^{}\right)\right\},`$
$`\sigma _q^{Q_k}(𝐫)`$ $`=`$ $`<[\widehat{\sigma }_q,\widehat{Q}_k]>`$
$`=`$ $`{\displaystyle \underset{nq}{}}\left(\phi _n^{}\widehat{\sigma }\widehat{Q}_k\phi _n^{}\widehat{Q}_k\phi _n^{}\widehat{\sigma }\phi _n^{}\right).`$
FIGURE CAPTIONS

Figure 1: Isovector E1 and isoscalar E2 giant resonances in $`{}_{}{}^{40}\mathrm{Ca}`$ and $`{}_{}{}^{208}\mathrm{Pb}`$ calculated with SkM\* forces in full RPA and SRPA at three levels of expansion as indicated. The strengths are weighted by Lorentz function with the averaging parameter $`\mathrm{\Gamma }=1.0`$ MeV.

Figure 2: Isovector E1 resonance in $`{}_{}{}^{40}\mathrm{Ca}`$ and $`{}_{}{}^{208}\mathrm{Pb}`$ calculated with SkM\* forces in full RPA and SRPA with the complete set $`k=14`$. The responses are weighted by Lorentz function with the small averaging parameter $`\mathrm{\Gamma }=0.1`$ MeV.

Figure 3: Isovector E1 resonance in $`{}_{}{}^{208}\mathrm{Pb}`$ calculated with SkI3 forces in full RPA and SRPA with a set of four operators. The responses are weighted by Lorentz function with averaging parameter $`\mathrm{\Gamma }=1.0`$ MeV.

Figure 4: Isovector E1 resonance in $`{}_{}{}^{208}\mathrm{Pb}`$ calculated in SRPA with the standard set of four operators and using the SkM\* force. Two different phases spaces are compared: full line = unrestricted space of excitations, dashed line = space restricted to $`ph`$ energies $`<30\mathrm{MeV}`$ as in previous figures. The responses are weighted by Lorentz function with averaging parameter $`\mathrm{\Gamma }=1.0`$ MeV.


format: text

text:
# Persistence of molecular excitations in metallic fullerides and their role in a possible metal to insulator transition at high temperatures

## I Introduction

The properties of strongly correlated materials have been in the focus of solid-state research for many years. To describe the competition between Coulomb and kinetic energies, the one band Hubbard model is widely used. However, many real systems exhibit orbital degeneracy, which adds a degree of freedom that is not always taken into account. Fullerides are one such example as the lowest unoccupied molecular orbital that forms the narrow t1u conduction band in the solid is triply degenerate. Together with the strong electron-phonon coupling characteristic of these materials, this leads to the possibility of Jahn-Teller distortions (JTD). Whereas predictions for the JTD can be done quite accurately for a single molecule Manini , their role in the solid is less clear. The broadening of the molecular levels should at first sight reduce the gain of energy associated to JTD. As it is estimated to be of the same order of magnitude than the gap opened by the JTD, the survival of these distortions in the solid can be questioned.

Yet, JTD appear to play a crucial role in many different fullerides. This paper concludes a serie of three papers devoted to NMR studies of different stoichiometries of alkali fullerides, where we have already seen evidence for such effects. It is widely believed that JTD contribute to turn Na2C60 and A4C60, which should be metals in a band picture, into non-magnetic insulators. In our first paper (called hereafter paper I BrouetPartI ), we have presented NMR data supporting this scenario and we refer the reader to references therein to the various papers invoking JTD in the properties of A4C60. Perhaps more surprisingly, they offer the most likely way for explaining the unexpected properties of the cubic quenched (CQ) phase of CsC60, the only alkali cubic metallic fulleride phase known so far besides A3C$`_{60}.`$ In this phase, we have shown in ref. BrouetPRL99  and paper II BrouetPartII  that the electronic properties are inhomogeneous on the local scale because spin-singlets are trapped on about 10% of the C60 balls. We believe that these singlets are stabilized by a JTD of the C60 ball which is energetically more favorable for C$`{}_{}{}^{2}{}_{60}{}^{}`$ than C$`{}_{}{}^{}{}_{60}{}^{}`$.

More generally, a simplified view of the effect of JTD on the t1u band as a function of its filling is sketched on Fig. 1. The gain of energy obtained from the JTD is found to be always larger for evenly charged C$`_{60}.`$ Then, JTD could induce attractive interactions for odd stoichiometries in order to promote the formation of the more stable C$`{}_{}{}^{2n}{}_{60}{}^{}`$ (hence, the presence of spin-singlet in CQ CsC60), whereas for even stoichiometries, they would induce repulsive interactions and favor localization (hence, the insulating Na2C60 and A4C60). This suggests an elegant way for rationalizing the different properties of even and odd stoichiometries Manini ; Heritier ; Gunnarsson  and the purpose of this paper is to determine to what extent this framework could be relevant for A3C60, extending ideas already presented in for the case of Na2CsC60 BrouetPRL2001 .

Despite the possible survival of these unusual molecular properties in the solid, there is a general consensus that the electronic properties of A3C60 can be understood in a rather conventional way, both for the metallic and superconducting states GunnarssonSupra . In the first part of this paper, we show that there are however in A3C60 deviations with respect to the conventional metallic NMR behavior, that we attribute to an enhanced stability of C$`{}_{}{}^{2}{}_{60}{}^{}`$ and C$`{}_{}{}^{4}{}_{60}{}^{}`$. The consequences of the formation of such “pairs” on the properties of the metal is an interesting issue to clarify. For this study, we choose two examples of superconducting fullerides Rb3C60 (T$`{}_{c}{}^{}=29`$ K) and Na2CsC60 (T$`{}_{c}{}^{}=12`$ K). While the structure of Rb3C60 is face centered cubic ($`fcc`$, space group Fm$`\overline{3}`$m) through the whole temperature range StephensNature91 , Na2CsC60 undergoes an orientational transition from $`fcc`$ to simple cubic ($`sc`$, space group Pa$`\overline{3})`$ around 300 K PrassidesScience94 . In the second part of this paper, we study the high temperature properties of two metallic fullerides where C$`{}_{}{}^{2}{}_{60}{}^{}`$ have been detected, Na2CsC60 and CsC60. The behavior departs from the one of a metal and can be best described by a progressive localization of the charge carriers with increasing temperature. This is likely related to a change of the role of the JTD as a function of temperature. However, the most important conclusion of this study might be that it reveals an overall similar behavior in these two compounds, which simplifies greatly the understanding of fullerides by unifying apparently conflicting properties. Up to now, CsC60 was thought to be metallic in its low temperature phase KosakaPRB95  and insulating in its high temperature phase TyckoPRB93 . The connection between these two behaviors was not understood, which has forbidden so far a comparison with A3C60 systems.

## II Detection of molecular excitations through spin-lattice relaxation measurements in superconducting fullerides

In this part, we demonstrate that the deviation of 1/T1T from the simple metallic behavior can be convincingly attributed to singlet-triplet distortions of JTD C$`{}_{}{}^{2}{}_{60}{}^{}`$ and C$`{}_{}{}^{4}{}_{60}{}^{}`$. To establish this, we discard in the first part other possible origin for this deviation. In the second part, we give a quantitative analysis of the coexistence between metallic relaxation and singlet-triplet excitations, which sets this model on a firm basis.

### II.1 Origin of the non-Korringa relaxation

In a metal, the relaxation of nuclear spins is usually dominated by their coupling with conduction electrons. This leads to a simple dependence of the relaxation rate on the density of states n(E$`{}_{f}{}^{})`$, known as the Korringa law Slichter .

$$\frac{1}{T_1T}=\frac{\pi k_B}{\mathrm{}}A^2n(E_f)^2$$
(1)
where $`A`$ is the hyperfine coupling (in erg) between 13C and conduction electrons. One anomaly in the NMR behavior of A3C60 is that 1/T1T frequently deviates from the 1/T1T = cst law that one consequently expects. An increase of 13C NMR 1/T1T with increasing temperatures has been reported for nearly all A3C60 compounds TyckoPRL92 ; ManiwaJapan93 ; ManiwaPRB95 ; StengerPRL95 , sometimes with moderate values (30 % in Rb3C60 TyckoPRL92 ) or very large ones (200% in Na2CsC60 ManiwaPRB95 ). To probe these aspects further, we have taken data in these two extreme cases in identical experimental conditions, with a 7 T applied magnetic field and standard saturation recovery pulse sequences, and over an extended temperature range. The results are presented on Fig. 2. In both compounds, the relaxation curves for the NMR magnetization are not exponential below 150 K, which causes a priori a problem for defining T1 at low temperatures. However, we have shown in paper I that the non-exponentiality is small and does not change between 10 and 150 K so that it cannot affect significantly these results. Moreover, the region below 150 K is precisely the one where the Korringa law is better obeyed. For the data presented here, T1 has been defined below 150 K as the mean value of a double exponential fit.

The narrow peak observed at 180 K in Na2CsC60 is due to a coupling with slowing down molecular motions, as was also observed in the isostructural Na2C60 BrouetPartI . As this shows that not only conduction electrons can contribute to the relaxation in our case, one could wonder whether similar peaks could be present (although not well resolved) in other temperature ranges or in Rb3C60. If so, could they explain part of the increase of 1/T1T ? The large temperature range for the data presented here allow to discard such a possibility, because these peaks would be symmetric while the increase of 1/T1T is essentially step-like. We refer the reader to our paper I for a detailed discussion of the characteristics of these molecular motion peaks.

From Eq. 1, some temperature dependence of 1/T1T could be expected if n(Ef) is not strictly constant as expected for a standard Pauli susceptibility. This could indeed happen, if there are narrow features in the density of states near the Fermi level. Those could be progressively disclosed as thermal expansion increases the lattice constant with increasing temperature. Because they are molecular solids, bound by weak Van der Waals interactions, fullerides are actually very compressible materials and the temperature dependence of n(Ef) must be seriously considered. Some studies have concluded that it is sufficient to explain the observed increase of 1/T1T for certain compounds StengerPRL95 . An obvious test to determine if the increase of 1/T1T is due to such an effect is to plot 1/T1T together with $`\chi ^2`$, as we have done for Na2CsC60 in ref. BrouetPRL2001 . We have found that the temperature dependence of $`\chi ^2`$ is much too small to explain that of 1/T1T, so that the deviation must be due, at least in this case, to the presence of an additional relaxation channel.

We want to reinforce this conclusion here by taking advantage of the comparison between Na2CsC$`_{60}`$and Rb3C60, over the large temperature range of the present experiments. The accuracy of the data presented on Fig. 2 at low temperature establishes that the increase is not regular. While 1/T1T is quite remarkably constant below 100 K (Na2CsC$`{}_{60}{}^{})`$ or 150 K (Rb3C60), it then increases up to room temperature, where it decreases (Na2CsC60) or saturates (Rb3C60). The lattice contraction on the other hand is expected to follow a smooth temperature dependence defined by the compressibility of the materials. Furthermore, the compressibility of $`sc`$ and $`fcc`$ phases are known to be quite similar Samara , so that the slightly different structure in the two compounds cannot explain such a difference. The thermal compressibility measured in Rb3C60 ($`\kappa =d(\mathrm{ln}a)/dT=3.10^5`$ K-1 Zhou ) corresponds to an increase of the lattice parameter $`a`$ by 0.16 Å between 10 and 400 K, very close from the one measured in Na2CsC60 (0.14 Å PrassidesScience94 ). One could argue that even though the temperature dependence of $`a`$ is the same, n(Ef) might display quite different variations with $`a`$ in $`sc`$ or $`fcc`$ materials YldirimSScom95 . This idea has been proposed to explain the different variation of $`T_c`$ with $`a`$ in Na2AC60 ($`sc`$) or A3C60 ($`fcc`$) materials (with A = K, Rb, Cs) because it is generally assumed that Tc depends on n(Ef) in a straightforward way according to BCS theory. However, the variation of Tc with pressure was later found the same in $`fcc`$ and $`sc`$ phases BrownPRB99 , so that the different variation of Tc is not related to n(Ef) but to an alkali effect. The comparison between Na2CsC60 and Rb3C60 is then legitimate and presents the advantage to cover a large variation in density of states. An even more serious drawback with such an explanation is that our study evidences a much larger increase in Na2CsC60 than in Rb3C$`_{60},`$ although a direct measurement of the susceptibility below 300 K shows the opposite behavior RobertPRB98 .

### II.2 Singlet-triplet excitations in A3C60

In Na2CsC60, the origin of the increase of 1/T1T was suggested by the comparison with Na2C$`_{60},`$ which exhibits a very similar relaxation behavior at high T despite it is insulating (see BrouetPRL2001  or Fig. 7). In this latter compound, the relaxation is due to singlet-triplet transitions between different JTD of a C$`{}_{}{}^{2}{}_{60}{}^{}.`$ Theoretical calculations have shown that for an isolated C$`{}_{}{}^{2}{}_{60}{}^{}`$, a JTD with a triplet ground state lies above the singlet JTD represented in Fig. 1 by Ea=140 meV, so that triplet states can be thermally populated Manini . Therefore, we have suggested that singlet-triplet excitations persist in Na2CsC60 because C$`{}_{}{}^{2}{}_{60}{}^{}`$ and C$`{}_{}{}^{4}{}_{60}{}^{}`$ are formed there on very short time scales (about 10-14 sec), which creates an additional relaxation channel explaining the increase of 1/T1T at high temperature. We now investigate if such an explanation could describe the evolution of the relaxation behavior between Na2CsC60 and Rb3C60.

The data led us to suggest that the relaxation can be divided between two relaxation channels, a metallic Korringa-like channel and a molecular channel corresponding to localized singlet-triplet excitations. For this latter term, we assume that the imaginary part of the susceptibility, the quantity measured by 1/T1 at the nuclear Larmor frequency (80 MHz for 13C in this measurement), can be described by a Lorentzian with a width 1/$`\tau `$, where $`\tau `$ is a lifetime characterizing the electronic excitations. This is is a usual assumption for localized magnetic moments, which yields for small $`\omega `$, 1/T$`{}_{1}{}^{}\chi \tau `$, where $`\chi `$ is the static paramagnetic susceptibility of the localized moment White . This leads to the following expression.

$$\frac{1}{T_1T}=\frac{k_B}{\mathrm{}}A^2\text{ }\left(\pi n(E_f)^2+\frac{\chi _{_{ST}}}{\mu _B^2}\frac{\tau _{st}}{\mathrm{}}\right)$$
(2)
where $`\tau _{st}`$ is the characteristic lifetime of a triplet state and $`\chi _{_{ST}}`$ an activated susceptibility describing the population of triplet states. The separation between the two terms is quite arbitrary, as the same electrons participate to both terms. Such a phenomenological decomposition should nevertheless capture the essential points of the physics as long as $`\tau _{st}`$ is shorter than the lifetime $`\tau _p`$ of a C$`{}_{}{}^{2n}{}_{60}{}^{}`$ in the metal. Some correlation between the parameters describing the “localized” states (like $`\tau _{st}`$) and the extended carriers states (like n(Ef)) can be expected as a result of this situation.

A source of inspiration for this correlation is the evolution of Rb4C60 from an insulator to a metal with applied pressure Kerkoud . This transition is not sharp, but a linear metallic term appears with increasing pressure that coexists with the localized activated one, as replotted on Fig. 3. As the metallic term grows, the molecular one can still be described by the same model, but the value of the gap and $`\tau _{st}`$ have to be reduced. This legitimates the assumption of coexisting molecular and metallic excitations, which is also observed in Na2C60 BrouetPartI . We believe that this property is actually essential to understand the originality of fullerides, in which molecular features are retained because the large electronic correlations forces one electron to spend a “long” time on each molecule before being transferred to the next.

The behavior in Rb4C60 tells us that $`\tau _{st}`$ directly depends on n(Ef). This suggests that the triplet states are relaxed by conduction electrons. Assuming that the triplet state can be treated as a local impurity on the timescale of its existence, we can use relaxation laws observed for localized impurities in a metal, that are analogous to a Korringa-like process White .

$$\frac{1}{\tau _{st}}=\frac{\pi }{\mathrm{}}\text{ }k_BT\text{ }J^2\text{ }n(E_f)^2$$
(3)
where J would be the coupling between the triplet and conduction electrons.

To test this model, which reduces the number of free parameters in Eq. 2, we have tried to reproduce the pressure evolution of the two terms in Rb4C60. To modelize the susceptibility associated with singlet-triplet excitations, we choose

$$\chi _{st}=\frac{8\mu _B^2}{k_BT}\frac{exp(\mathrm{\Delta }/T)}{2+3\mathrm{exp}(\mathrm{\Delta }/T)}$$
(4)
which was found in reasonable agreement with the susceptibility of Na2C60 BrouetPartI . The coupling $`J`$ is treated as a shared parameter between the three sets of data and the best fits are found for $`J=10^{13}`$ erg. Fits using equations 2, 3 and 4 with this value for J are shown on Fig. 3. They are in fair agreement with the data, which supports this model. The density of states and the value of the gap used in each case are displayed in Table 1. This model can be applied as well to Na2C$`_{60}.`$ Fig. 3 shows that, with a larger gap value already noticed previously BrouetPRL2001  but keeping the same value for J, the data are also well reproduced.

Does this model also apply to the case of A3C60 ? Eq. 3 implies that the increase of 1/T1T due to the localized term should show up more clearly when n(Ef) is small, which is in qualitative agreement with the evolution of 1/T1T between Na2CsC60 and Rb3C60. Remarkably, the previous fitting procedure with the same J value gives a correct order of magnitude for the increase of 1/T1T in both compounds. In this model, the magnitude of the high temperature increase of 1/T1T is uniquely fixed by its value at low temperature. The best fits displayed on Fig. 2 with the parameters given in table 1 show an excellent agreement in the case of Rb3C60 but is poorer for Na2CsC$`_{60}.`$ In this latter compound, it seems that a different regime is relevant for high temperatures, which will be discussed in the next section. Even below 300 K, the increase of 1/T1T cannot be ascribed to a single gap value. Similar problems were encountered in the case of Na2C60, which led us to suggest in paper I that the equilibrium between different JTD could be modified by changes in details of the structure as a function of temperature.

To check the consistency of our analysis, the triplet relaxation time $`\tau _{st}`$ should be compared to the lifetime $`\tau _p`$ of the C$`{}_{}{}^{2n}{}_{60}{}^{}`$. This latter time can only be longer that the average time spent by one electron in the vicinity of a C60 ball $`\tau _m`$. A lower bound for $`\tau _m`$ can be estimated by the time required for an electron to travel at the Fermi velocity from one ball to the other, $`\tau _m`$= a / v$`{}_{f}{}^{}\mathrm{}n(E_f`$). This is calculated in Table 1 and the comparison with $`\tau _{st}`$ deduced from the fits shows that they have similar order of magnitude at room temperature, so that this analysis is consistent. We note that we have used here the susceptibility given by Eq. 4, although one could assume that only a fraction of it, corresponding to the actual concentration of C$`{}_{}{}^{2n}{}_{60}{}^{}`$ in A3C60, should be used. This would give slightly longer value for $`\tau _{st}`$.

Therefore, this simple model gives an overall satisfying description of the relaxation behavior in these compounds. The correlation between the molecular and metallic term might however be even more intricate. We have assumed that the JTD C$`{}_{}{}^{2}{}_{60}{}^{}`$ can be considered as an isolated molecular entity, even if it is restricted to a limited timescale. The fact that the gap value is not fixed, but tends to decrease with increasing n(Ef) in compounds with identical structures like Na2C60 and Na2CsC$`_{60},`$ shows explicitely that there is a further correlation between the two terms. $`\mathrm{\Delta }`$ is not purely a molecular value and a more sophisticated description, introducing for example a screening of the C$`{}_{}{}^{2n}{}_{60}{}^{}`$ by conduction electrons, would be required to describe completely this behavior.

In addition, we have mainly tried to explain so far the increase of 1/T1T which roughly takes place between 200 and 300 K. At higher temperatures, Eq. 3 predicts a decrease of the lifetime of the triplet states as 1/T. This could change the balance between the two relaxation channels as a function of temperature. The simple estimate that we have used for $`\tau _m`$ is also likely to break down at high temperature because it requires that the mean free path for the electronic motion $`l`$ is longer than $`a`$. On the contrary, it is known that very small values for the mean free path are deduced from resistivity measurements for fullerides at high temperature, like 1-2 Å at 500 K in Rb3C$`_{60},`$ HebardPRB93 . Therefore, we could expect complications at high temperature, that might reveal, on the other hand, valuable information on the charge transport in these materials and we turn to this topic in the next section.

## III Are “metallic fullerides” metallic up to high temperatures ?

Some anomalies in the high temperature behavior of fullerides have been noticed for a long time. In Rb3C$`_{60},`$ the photoemission spectra at high temperature does not display a clear Fermi edge Golden . In Na2CsC60, optical conductivity measurements indicate a disappearance of the Drude-like peak above 300 K Cegar . All this motivates us to investigate the high temperature region in more details and we indeed report here two cases, Na2CsC60 and CsC60, where the high temperature behaviors are far from that of a simple metal and are even suggestive of a localization of the charge carriers. For all the data reported here, we have checked that the behavior is fully reversible with temperature, ensuring that there is no deterioration of the sample quality and/or stoichiometrie.

### III.1 Na2CsC60

#### III.1.1 progressive localization of the charge carriers inferred from 1/T1 behavior

Fig. 2 shows that 1/T1T decreases regularly from 300 to 700 K in the $`fcc`$ phase of Na2CsC60. This is anomalous as among the two terms identified in the relaxation so far, the metallic and the singlet-triplet components, the first one is expected to yield a constant value and the second one an increase of 1/T1T or a slight decrease depending on the precise value of the gap (see the fit of Fig. 2). In Rb3C$`_{60},`$ there might be a similar anomaly above 500 K, where 1/T1T is somewhat lower than the fitted curve, but this cannot be concluded unambiguously within experimental accuracy and in the following we focus on Na2CsC$`_{60}.`$

As the use of 1/T1T might introduce an artificial bias by emphasizing the metallic component, we have replotted on Fig. 4, 1/T1 for 13C NMR, together with the result from 133Cs NMR, which displays essentially the same temperature dependence, except for the molecular motion peak. In our case, the interpretation of these data is not straightforward as we do not know how to separate the metallic term from the ST component. Within the metallic framework of the Korringa law, a decrease of 1/T1T would be assigned to a reduction in the density of states signaling the onset of a broad metal to insulator transition. If it is a change in the singlet-triplet component, it could be due to modifications of the apparent value of the gap as was observed in Na2C$`_{60}.`$ To find a way out of this problem, a comparison with the ESR susceptibility shown in the inset of Fig. 4 is helpful. In the susceptibility, both the metallic and singlet-triplet terms are also present but with different ratios as they are not weighted by the lifetime of the excitations, like 1/T1 is. $`\chi `$ is dominated at low T by a constant Pauli component and then increases up to 350 K, in our opinion mainly as a result of singlet-triplet excitations. Like 1/T1T, it is clearly anomalous above room temperature, where it starts decreasing sharply. This is very different from the behavior of $`\chi `$ in Na2C60, where it keeps increasing at high temperature, which tells us that a modification of the ST component alone cannot account for this result. Above 350 K, we observe a strong divergence in behavior between $`\chi `$ and 1/T$`_1,`$ $`\chi `$ falls far below the metallic value, whereas 1/T1T is still far above. This can only be understood if the lifetime of the excitations diverges above room temperature. This reminds of a situation observed at the transition between solid and liquid in X2Te3 (X=In, Ga, Sb) Warren . Upon melting, the increasing disorder induces a progressive localization of the charge carriers that leads to a decrease of the Pauli susceptibility, due to the reduction of the carrier density. At the same time, 1/T1T is found to increase, as the lifetimes of the excitations increases rapidly. To account for this effect n(Ef)2 in Eq. 2 has been replaced by the authors of ref Warren , by n(Ef)$`\tau _m/\mathrm{},`$ where $`\tau _m`$ is now free to deviate from the metallic value $`\mathrm{}n(E_f`$). In our case, both $`\tau _m`$ and/or $`\tau _{st}`$ could be responsible for the increase of 1/T1T. However, in both cases, we can conclude that there is a progressive localization of charge carriers.

#### III.1.2 nature of the localized charge carriers inferred from the static NMR alkali spectra

If the localization process went as far as to correspond to a static charge separation, we could expect to observe changes in the NMR spectra as a result of inequivalent alkali sites neighboring a different number of C$`{}_{}{}^{2n}{}_{60}{}^{}`$. This is the case in CQ CsC60 where a splitting of the Cs spectrum is observed at low temperature, corresponding to the localization of 10% C$`{}_{}{}^{2}{}_{60}{}^{}`$, as discussed in details in paper II. In Na2CsC60, static charge separation can be ruled out because 133Cs and 23Na spectra consists of one featureless narrow line (with the exception of the Na “T’ ” line). This is consistent with the value given previously for $`\tau _{st}`$ which is much shorter than the NMR time scale (a few ms corresponding to the inverse width of the spectrum), so that a motional narrowing of the spectra should still take place even if $`\tau _{st}`$ increases by many order of magnitudes.

There are however clear anomalies in the alkali NMR which are revealed when comparing the shifts $`K`$ for 133Cs and 23Na lines, presented on Fig. 5, together with the ESR susceptibility. A simple scaling between $`K`$ and $`\chi `$ is expected, $`K=\sigma +A\chi `$, where A is the hyperfine coupling and $`\sigma `$ a reference chemical shift. On Fig. 5, it can be seen that both nuclei sense an increase of the susceptibility between 150 K and 300 K, but they do not reproduce similarly the decrease seen by ESR above 300 K. This also suggests, although in a puzzling way, that the properties of the material are changing above room temperature. Let us now discuss quantitatively the behavior of the shift for the two nuclei.

##### 23Na NMR

In the $`sc`$ phase, we can scale the 23Na shift and the susceptibility, as shown by the line of Fig. 5a by using $`\sigma =300`$ ppm and A = 4500 Oe/$`\mu _B`$. It is interesting to compare these numbers with those found for 23Na NMR in paper I for the isostructural Na2C$`_{60},`$ $`\sigma `$=-65 ppm and $`A`$= 2300 Oe/$`\mu _B`$. The difference for $`\sigma `$ is quite large and not expected as $`\sigma `$ essentially depends on the diamagnetism of the Na+ filled orbitals which should not change much. To avoid this problem, it might be more realistic to assume that the hyperfine coupling to the metallic component is smaller than the one to localized excitations, maybe because Na moves a little bit towards the C$`{}_{}{}^{2}{}_{60}{}^{}`$ or C$`{}_{}{}^{4}{}_{60}{}^{},`$ which increases the hyperfine coupling for this term. This automatically gives a value of $`\sigma `$ closer to that of Na2C60, as the shift is dominated at low T by the metallic term.

Through the $`sc`$-$`fcc`$ transition at 310 K, we have observed in Na2C60 an increase by 30% of the hyperfine coupling probably associated to the change of structural environment of Na (see paper I). This is also clearly present in Na2CsC60. As a guide to the eye, we have plotted the variation of $`K`$ that would be expected for 23Na NMR from the one of $`\chi `$ assuming no change for $`\sigma `$ and a continuous value for $`\chi .`$ This only emphasizes what was immediately clear, $`K`$ does not follow the static susceptibility at high temperature, it is nearly constant while $`\chi `$ drops by a factor 2. Relying on the previous findings i) the characteristic lifetime of the localized excitations starts to increase above 300 K and ii) the hyperfine coupling to the localized excitations is particularly strong, we would like to conclude that the reason for the different behavior of K and $`\chi `$ at high temperature is due to an increase of the hyperfine coupling associated with the progressive localization of the carriers. Indeed, we know from CQ-CsC$`_{60},`$ studied in paper II, that the coupling can be very strong with localized singlets in a metallic environment. This enhancement disappears if all the balls are equivalent and bear a singlet, like in Na2C60, because Na has no reason to move towards a particular ball.

##### 133Cs NMR

The 133Cs shift presents a similar behavior but it displays a significant temperature variation below 100 K, which does not appear in the ESR data. Such a linear temperature dependence of the shift has already been observed in the NMR of various alkali in fullerides AlloulPhysica ; Stenger ; ZimmerEurophysics . Its origin is not well understood, but as it is not related to the electronic susceptibility AlloulPhysica , it is probably of orbital origin and we will include this linear variation in the chemical shift $`\sigma `$. Choosing $`\sigma =200+0.25`$ $`T,`$ as indicated on Fig. 5 by the dotted line, allows to extract an additional increase, clearly present between 150K and 300K which scales with the susceptibility with A= 5000 Oe/$`\mu _B`$. As done for 23Na, we then extrapolate the variation of $`K`$ at high temperatures assuming that $`\sigma `$ and $`\chi `$ do not change, which requires here a reduction of A in the $`fcc`$ phase to 4500 Oe/$`\mu _B`$ rq1 . Hereagain, there is a significative difference between K and $`\chi ,`$ although less dramatic than for Na. We also observe that the scaling between K and $`\chi `$ is not lost at the structural transition but at somewhat higher temperatures T$``$ 400 K, making it unlikely that the structural transition could be responsible for the change of the electronic properties of Na2CsC$`_{60}.`$ If the divergence between K and $`\chi `$ is to be attributed again to an increase of the hyperfine coupling, this yields A=6000 Oe/$`\mu _B`$ at 650K, which corresponds to a 25 % increase from 350 K to 650 K. Such an increase should also enhance 1/T$`_1,`$ proportional to A$`^2,`$ by a factor 1.5. On Fig. 4, we see that 133Cs 1/T1 indeed increases a little bit more rapidly than 13C 1/T1 at high temperatures. As 13C cannot move, its hyperfine coupling can be used as a reference for the “unenhanced” value of 1/T$`_1.`$ This predicts 1/T1=13 sec-1 instead of the observed 1/T1=17 sec$`^1,`$ which is close from a factor 1.5 and supports the idea of an enhancement of the hyperfine coupling rq2 .

In summary, many anomalies appear above 400 K in Na2CsC60, which gives evidence for a change in the electronic properties. Although, the data are conflicting at first sight, because different probes measure apparently different behaviors, we propose that this can be understood assuming that we enter an intermediate regime between metallic and localized behavior. In this regime, 1/T1 is dominated by the increase of the lifetime of the charge carriers around a particular C60 ball, which also modulates the value of the alkali hyperfine couplings. We cannot decide unambiguously whether the charge carriers are localized as C$`{}_{}{}^{3}{}_{60}{}^{}`$ or C$`{}_{}{}^{2n}{}_{60}{}^{}`$. Both are consistent with the 1/T1 measurements. The change in the alkali hyperfine coupling somewhat favors C$`{}_{}{}^{2n}{}_{60}{}^{}`$, as we could expect a larger motion of the alkali towards the differently charged C60. The idea that a motion of the alkali ions, especially of the light Na ion, is playing an important role in the high temperature properties of this phase was already suggested in Cegar . However, we cannot rule out that the appearance of a localized charge is sufficient to trigger such a motion and modify the alkali hyperfine coupling regardless of the charge value.

Because we have detected in the previous part of this paper an important role of JTD in the metallic phase, it seems reasonable to associate the progressive localization at high temperatures with a conflict between the stabilization of JTD and the delocalization of the charge carriers. We can for example speculate that the population of many different distortions at high temperatures introduces disorder for the electronic system, which hinders the electronic motion. These electrons could then be localized as C$`{}_{}{}^{3}{}_{60}{}^{}`$ or C$`{}_{}{}^{2n}{}_{60}{}^{}`$ depending on the precise balance between Coulomb repulsion and Jahn-Teller attraction, but, more importantly, this behavior reveals by defect that a cooperation is required between JTD and electronic motion to obtain the metallic state. This agrees with the idea that JTD, at last at low temperatures, help the systems with odd stoichiometries to become metallic by stabilizing the C$`{}_{}{}^{2n}{}_{60}{}^{}`$, which have to be formed as a byproduct of electronic jumps.

### III.2 The high temperature cubic phase of CsC60

Another very anomalous high-temperature phase is the one of CsC60. It appears to have a Curie-like susceptibility Chauvet , the 133Cs NMR shift also follows a Curie law within experimental accuracy and 1/T1 is nearly constant as expected for a paramagnetic insulator TyckoPRB93 . All of this was taken as evidence for electronic localization. As for NMR, we have checked that this behavior holds up to 700 K for the 133Cs shift (see Fig. 6) and 13C 1/T1 (see Fig. 7, actually T$`{}_{}{}^{1}{}_{1}{}^{}`$ is not really constant but decreases by 25 % from 400 K to 700 K).

Below 350 K, a structural transition takes place to a polymerized phase with very different properties. However, the cubic phase can be studied at low temperatures (below 130 K) by quenching the high temperature phase and it seems to be metallic. A Pauli susceptibility KosakaPRB95  and a Korringa law for 13C NMR above 50 K BrouetPRL99  are observed. The behavior of the susceptibility in the two phases are reported on Fig. 6 to illustrate this point. The difference is very puzzling, as the structure is the same except for the orientational order (the CQ phase of CsC60 is $`sc`$ Lappas , while the high temperature phase is $`fcc`$). If this structural transition were to play a role in the conductivity, the high temperature phase, which is structurally equivalent to metallic A3C$`_{60},`$ should be the more directly comparable to other fullerides, whereas it is this phase that displays “anomalous” properties. This distinct behavior in the two phases casts doubt on whether 1 electron in the t1u band yields a metallic or insulating behavior.

On the other hand, it appears now quite similar to that of Na2CsC60 and it seems less difficult to connect the properties of the two CsC60 phases. Na2CsC60 also resembles a good metal at low temperature and, for some properties, particularly the ESR “Curie-like” susceptibility, an insulator above 300 K. If the data between 100 and 300 K were missing for Na2CsC60, one would also conclude that we are dealing with two entirely different phases. We have seen through this paper that including singlet-triplet excitations of JTD balls offers a convincing way to explain this apparent discrepancy, and the same could be true for CsC$`_{60}.`$ The underlying similarity is perhaps best illustrated by looking at the temperature dependence of 1/T1 for 13C in various phases, as we do on Fig. 7. In this context, it is tempting to interpret the gap between the values of 1/T1 in the CQ phase and those of the high temperature phase by an increase due to a singlet-triplet component. In CQ CsC60, just before the transition to the dimer phase, we do observe an increase of 1/T$`_1,`$ which might correspond to such a component. Adjusting a fit to Eq. 2 to these points and the high temperature value gives the line sketched on Fig. 7. This yields a smaller gap $`\mathrm{\Delta }=50`$ meV than for the other compounds, but this is consistent with the fact that 1/T1 decreases at high temperatures. Indeed, Fig. 7 shows a systematic relation between the value of the gap and the behavior at high temperature: the lower the gap, the higher 1/T1 and the smaller (eventually negative) the slope of its variation at high temperature. Qualitatively, this is not surprising as, with a small gap, the susceptibility is dominated at high temperature by the Curie law for the triplet states (see Eq. 4).

Putting these data together then gives a fairly good understanding of the CsC60 system. We know from our study of CQ CsC60 that JTD C$`{}_{}{}^{2}{}_{60}{}^{}`$ are formed at low temperatures but, presumably because there is no symmetric JTD for neutral C$`_{60},`$ they get localized on a small fraction of the C60 balls, contrary to C$`{}_{}{}^{2n}{}_{60}{}^{}`$ in superconductors. The singlets start to move as the temperature increases and the lifetime of a C$`{}_{}{}^{2}{}_{60}{}^{}`$ in the metal decreases exponentially as shown in paper II. We can speculate that, if it were stable, the properties of cubic CsC60 would be very similar to those of a metal like Na2CsC60 between 100 K and 300 K and would evolve gradually towards an almost insulating state at high temperatures, which is the one observed above 400 K.

## IV Conclusion

To conclude, we have proposed in this paper that the metallic character of fullerides with an odd number of electrons per C60 is driven by the formation of pairs of electrons on very short time scales stabilized by Jahn-Teller distortions. We argue that the deviation of the NMR 1/T1T from the canonical “metallic” Korringa law observed in many A3C60 systems can be convincingly attributed to the presence of singlet-triplet excitations of these electronic pairs. By studying two superconducting metals with very different density of states, we evidence a dependence of the characteristic time $`\tau _{st}`$ of the triplet state for these pairs on n(Ef) that makes the detection of singlet-triplet excitations easier in compounds with small density of states. This explains why the effect has not be clearly identified up to now, the molecular excitations being masked in the most heavily studied compounds like K3C60 and Rb3C60.

We have then presented new data at high temperatures that show a clear evolution of the behavior of the electronic properties in at least two compounds Na2CsC60 and CsC60. In the first case, we interpret the breakdown of the scaling between various quantities at about 400 K as the onset of a progressive localization of the charge carriers. We suggest that this localization is due to a change of the status of JTD in the metal as the temperature increases. The localized charge carriers could either be the electronic pairs formed via JTD or independent electrons, depending on the precise nature of the high temperature evolution. A similar, probably slightly more efficient, localization takes place in CsC60. Although this phase was already believed to be insulating, seeing this as a metal to insulator transition allow to connect its properties to those of the low temperature metallic phase. This adds greatly to the understanding of these systems, as all metallic phases can now be described within the same framework. JTD appear as a key ingredient in this behavior, which makes it likely that the difference between even and odd stoichiometries, respectively insulating and metallic, is due to the sensitivity of Jahn-Teller distortions to the parity of the C60 charge. As recalled in Fig. 1, JTD are always a little bit more stable for even parity.

This work also raise new questions on the nature of charge transport in A3C$`_{60}.`$ If the enhanced lifetime of C$`{}_{}{}^{2}{}_{60}{}^{}`$ and C$`{}_{}{}^{4}{}_{60}{}^{}`$ is a general feature of A3C$`_{60},`$ it is of course important to examine their role within the metal. The role that we assigned here to Jahn-Teller distortion is an example of how electronic properties are coupled to phonons and could be a first step towards the concept of polarons. However, we do not know how fast the JTD can adjust itself to a change of the charge of the molecule, produced by jumps of electrons from one ball to the other. The idea of a polaron would imply that one electron moves with a given JTD. In our case, we could as well assume that a JTD is fixed for a given molecule, regardless of its charge, which only increases the lifetime of favorable JTD C$`{}_{}{}^{n}{}_{60}{}^{}`$ configurations (namely n even). The lack of experimental knowledge about the time scale for the Jahn-Teller distortion, relatively to the motion of electrons, makes it difficult to discriminate between these possibilities. Nevertheless, some correlation between them is required to establish a coherent band-like charge transport. This could be destroyed with increasing temperature resulting in hopping like process. This might be what we witness in Na2CsC$`_{60}`$ and this bears interesting similarities to the metal to insulator transition observed at high temperature in manganites with ferromagnetic metallic ground states. In these systems also strong electronic correlations are coexisting with the possibility of JTD of the doubly degenerate eg orbitals. It has been proposed that the occurrence of JTD at high temperatures, which forbids hopping in this case, adds to the double exchange mechanism to trigger a transition to the high temperature insulating phase Millis .