text:
\section{INTRODUCTION}

Few--body nuclei with a number of nucleons A between 5 and 16
are a particularly interesting testground for nuclear theory.
They lie in the range between the classical few--nucleon systems ($A \le 4$)
and the smallest nuclei, which can be described realistically
starting from a mean field ansatz. Therefore one hopes that these nuclei could
build a link between few--body and many--body physics. Presently quite an
effort is made for a better understanding of these intermediate systems.
Specific interest is devoted to halo--nuclei, but also  the less exotic
nuclei of this mass range are investigated thoroughly. Many theoretical
techniques for the calculation of their ground states have been
imported from the classical few--body field, where there has been a
considerable
progress in the last years. In fact for the classical few--body systems
rather different approaches have been developed and proven to lead to precise
results. These methods include solutions of the Faddeev--Yakubovsky
equation, variational (VMC) and Green Function
Monte Carlo (GFMC), the Hyperspherical Harmonic (HH) ansatz,
the stochastic variational  (SVM), as well as coupled cluster and
resonating group methods.

For nuclei with $A>4$ a similar level of precision has not yet been
reached. Here GFMC has led to the most accurate results.
Exact bound state energies with realistic $NN$ interactions have been
calculated for $A \le 9$ ~\cite{PPCPW98}.
Unfortunately wave functions cannot be generated with this method. Recently,
a very powerful tool to calculate few--body wave functions has been developed
with the SVM ansatz ~\cite{{Varga95},{Varga96}}.
However, this approach is probably most suitable for
systems with $A < 8$. For nuclei with $A \ge 8$ rather good results have been
obtained in the Integro--Differential Equation Approach (IDEA)
~\cite{Brizzi96}, which uses the
HH expansion at first order, and with a variational method, the Translational
Invariant Configuration Interaction Method (TICI) ~\cite{{BFBBG90},{GMNBPW96}},
which is inspired by the coupled cluster method.

In this work we present the general formulation of the method, which combining
the main ideas of the HH expansion, the pair correlation ansatz and the
Faddeev approach allows to calculate wave functions of few--body systems.
We apply it to calculate the binding energies of
nuclear systems with 6, 8 and 12 particles.
Including higher order HH functions we make the
nontrivial step beyond the IDEA approach.
The difficulty in constructing antisymmetric $A$--particle HH is overcome
by the use of a recently developed very efficient recursive method
~\cite{Nir97}, where HH basis functions are constructed,
belonging to well defined
irreducible representations of the orthogonal and symmetric groups.

The paper is organized as follows. The method is described in Section
\ref{gen}, while the construction of the HH basis function is
briefly reviewed in Section \ref{basis}. Section \ref{matrix} illustrates
how the matrix elements are calculated and the obtained results for the
binding energies of nuclei with $A = 6, 8$ and 12 are discussed in
Section \ref{results}. Conclusions are drawn in Section \ref{concl}.

\section{GENERAL FORMULATION OF THE METHOD}\label{gen}

Our aim is the solution of the Schr\"odinger equation for a system of
A particles interacting via a two--body potential. After subtraction of the
center of mass Hamiltonian the problem reduces to the solution of the
internal Hamiltonian in terms of $N=(A-1)$ Jacobi coordinates $\vec\eta_i$.
One can  write the wave function in terms of an  HH expansion
\begin{eqnarray} \label{HH}
    \Phi(\vec\eta_1..\vec\eta_{A-1}, s_1..s_A, t_1..t_A)&=&\sum_{K\,\nu}
              R_{K\,\nu}(\rho) H_{K\,\nu}(\Omega,s_1..s_A, t_1..t_A)\cr
    &\equiv& \sum_{K\nu}\Phi_{K\nu}(\vec\eta_i,s_1..s_A, t_1..t_A)\,,
\end{eqnarray}
where, in case of nucleons,
\begin{eqnarray} \label{H}
  H_{K\,\nu}(\Omega,s_1..s_A, t_1..t_A) =
         \sum_{Y_{A-1}}
         \frac{\Lambda_{\Gamma_{A},Y_{A-1}}}{\sqrt{| \Gamma_{A}|}} \,
         &{\cal Y}&_{K_N L_N M_N
          \, \Gamma_{A} Y_{A-1}\, \alpha^K_N}(\Omega)\,\times\cr
          &{\cal X}&_{S S^z T T^z \, \tilde{\Gamma}_{A},\tilde Y_{A-1}
         \, \alpha^{ST}_A}(s_1..s_A, t_1..t_A)
\end{eqnarray}
represent the hyperspherical harmonic functions coupled to the spin--isospin
basis functions to yield the totally antisymmetric wave function.
The ${\cal Y}_{K_N L_N M_N \, \Gamma_{A}Y_{A-1} \, \alpha^K_N}(\Omega)$
are the so called symmetrized hyperspherical functions that depend
on the $N$ Jacobi coordinates.
They are
HH functions with hyperspherical angular momentum $K=K_N$, and good angular
momentum quantum numbers $L_N, M_N$ that belong
to well defined irreducible representations (irreps)
$\Gamma_{1} \in \Gamma_2 \ldots \in \Gamma_A $, of the permutation
group--subgroup chain
${\cal S}_1 \subset {\cal S}_2 \ldots \subset {\cal S}_A $,
denoted by the Yamanouchi symbol
$[ \Gamma_A, Y_{A-1} ] \equiv [ \Gamma_A,\Gamma_{A-1},\ldots,\Gamma_1 ]$.
The dimension of the irrep $\Gamma_{n}$ is denoted by $| \Gamma_{n}|$
and $\Lambda_{\Gamma_{A},Y_{A-1}}$ is a phase factor \cite{Akiva88}.
Similarly, the function
${\cal X}_{S S^z T T^z \, \tilde{\Gamma}_{A}\tilde{Y}_{A-1}
         \, \alpha^{ST}_A}$
is a symmetrized spin--isospin state with good quantum numbers $S, S^z, T$
and $T^z$.
The label $\alpha^K_N$ ($\alpha^{ST}_A$) is needed to remove the degeneracy
of the hyperspherical (spin--isospin) states with a given symmetry.
The argument of the hyperradial function  $R_{K\,\nu}$
is the hyperradius $\rho$ defined by
\begin{equation}
\rho^2=\sum_i \eta_i^2\,,
\end{equation}
$\Omega$ is the (3A-4)--dimensional hyperangle, and for brevity we shall
use $\nu$ for all the quantum numbers but $K$, i.e.
$\nu \equiv (L_N M_N \, S S^z T T^z \, \Gamma_A \alpha^K_N \alpha^{ST}_A)$.
In what follows we shall use the subscript $N$ to denote $A$ particle
quantum numbers that depend on the Jacobi coordinates, and the subscript
$A-1$ for the $A-1$ particle subsystem.

It is well known that finding a solution of the Schr\"odinger equation in
terms of the ansatz (\ref{HH})
can be very difficult because the number of basis functions increases very
fast with $K$ and in order to
have a real convergence one must use a huge number of basis functions
\cite{Kievsky97}.
Therefore a correlation function is advantageous to give the wave function
a proper behavior ~\cite{FE72}.
Its advantages have been extensively verified for classical few-body
systems reaching a high level of accuracy \cite{Kievsky93,Victor89}.
A general ansatz within the two--body correlation scheme is the Jastrow
factor
\begin{equation} \label{Jastrow}
   \Psi = \prod_{i<j} f_{ij} \Phi\,,
\end{equation}
where $f_{ij}$ is a two--body correlation function.
However, the use of the Jastrow ansatz leads to $3A-3$
dimensional integrals. Therefore it is more convenient to use the so called
pair correlation ansatz
\begin{equation} \label{Psi}
\Psi=\sum_{i<j} \chi_{ij} \Phi,
\end{equation}
because in this case one can use the Faddeev approach which leads to at
most four--body integrals.
In the Faddeev approach the Schr\"odinger equation is replaced by
equivalent equations,
\begin{equation} \label{Faddeev}
   (T-E)\Psi_{ij}=-V_{ij} \Psi \; ,
\end{equation}
where $ \Psi = \sum_{i<j} \Psi_{ij}$.
In order to speed up the convergence these equations can be further
modified \cite{Brizzi96} to include
the contribution of the hypercentral potential explicitly.
The hypercentral potential $V_{hc}(\rho)$ is defined as the projection
of the two--body
interaction on the subspace of the lowest order hyperspherical state, i.e.
the hyperspherical state expressed in  Eq. (\ref{H}), with the minimal
$K$:
\begin{equation} \label{Vhc}
  V_{hc}(\rho)=
        \int d \Omega \, H_{K_{min} \nu}^{\dagger}(\Omega,s_i,t_i)
         V(\vec{r}_{1,2},s_1,s_2,t_1,t_2)
         H_{K_{min}\nu} (\Omega,s_j,t_j)\,.
\end{equation}
Here the integration is carried over the $3N-1$ dimensional hypersphere,
and  an implicit summation over all the spin--isospin states is understood.
With the help of this definition we can rewrite Eq. (\ref{Faddeev}) as
follows
\begin{equation} \label{Fad_hc}
   \left[ T+\frac{A(A-1)}{2}V_{hc}(\rho)-E \right]\Psi_{ij}=
   -[V_{ij}-V_{hc}(\rho)] \Psi \; .
\end{equation}
Motivated by the pair--correlation ansatz (\ref{Psi}) we shall expand
the Faddeev amplitude $\Psi_{ij}$ in the following way,
\begin{equation} \label{Psi_ij}
  \Psi_{ij}=\sum_{p} \Phi_p (\rho,\Omega) \chi_p(z_{ij})
	    =\sum_{K \nu p} R_{K\nu p}(\rho)
               H_{K\nu} (\Omega,s_1..s_A, t_1..t_A)\chi_p(z_{ij}) \,,
\end{equation}
where $\chi_p(z_{ij})$ is a polynomial of order $p$, and
$z_{ij}$ is related to the relative two--body distance through
\begin{equation}
	z_{ij} = r_{ij}/\rho \,.
\end{equation}
Although expanding the correlation function in terms of $z_{ij}$
should be equivalent to an expansion in $r_{ij}$, it is more stable
numerically to take $\chi_p(z_{ij})$ since $r_{ij}$ goes to zero with $\rho$.
Substituting the expansion of Eq. (\ref{Psi_ij}) into Eq. (\ref{Faddeev})
we get
\begin{equation} \label{Fad2}
   \left[-{1 \over 2}\sum_{n=1}^{A-1} \Delta_n
   +\frac{A(A-1)}{2}V_{hc}(\rho) -E \right]
   \sum_{p} \Phi_{p} \chi_p(z_{ij})
   =-[V_{ij}-V_{hc}(\rho)] \sum_{p} \Phi_{p} \sum_{k>l} \chi_p(z_{kl}) \,,
\end{equation}
where $\Delta_n$ is the Laplace operator associated with the $n$'th Jacobi
coordinate $\vec{\eta}_n$.
Multiplying  Eq. (\ref{Fad2}) on the left by
$( H_{K'\nu'} \chi_{p'}(z_{ij}) )^{\dagger}$ and using the notation
\begin{equation} \label{defME}
   _{ij}< K' \nu' p' | \hat O | K \nu p >_{kl} =
   \int d \Omega \, H_{K'\nu'}^{\dagger}(\Omega,s_n,t_n)
   \chi_{p'}^*(z_{ij})
   \hat O H_{K \nu}(\Omega,s_m,t_m) \chi_p(z_{kl})
\end{equation}
we get
\begin{eqnarray} \label{Fad3}
 & & \sum_{K \nu p} \,
    _{ij}< K' \nu' p' |-\frac{1}{2}\sum_{n=1}^{A-1} \Delta_n
    +\frac{A(A-1)}{2}V_{hc}(\rho) - E  | K \nu p >_{ij}
    R_{K \nu p}(\rho) = \nonumber \\ & & \hskip 2truecm
   -\sum_{K'' \nu'' p''} \, _{ij}< K' \nu' p' | [V_{ij}-V_{hc}(\rho)]
    \sum_{k>l}| K'' \nu'' p'' >_{kl} R_{K'' \nu'' p''}(\rho) \; .
\end{eqnarray}
If the norm matrix
\begin{equation}
  A_{K' \nu' p', K \nu p} \equiv
     _{ij}< K' \nu' p' | K \nu p >_{ij}
\end{equation}
is non singular, one can define the projection operator
\begin{equation}
  \hat Q = \sum_{K \nu p, K' \nu' p' } | K' \nu' p' >
           A_{K' \nu' p', K \nu p}^{-1} < K \nu p | \,.
\end{equation}
The projection operator commutes with $V_{ij}$ and
one can rewrite Eq. (\ref{Fad3}) in the following form
\begin{eqnarray} \label{Fad-A}
 & & \sum_{K \nu p} \,
    _{ij}< K' \nu' p' |-\frac{1}{2}\sum_{n=1}^{A-1} \Delta_n
    +\frac{A(A-1)}{2}V_{hc}(\rho)    - E  | K \nu p >_{ij}
    R_{K \nu p}(\rho) = \nonumber \\ & & \hskip 2truecm -
    \sum_{K \nu p} \sum_{K'' \nu'' p''} \sum_{K''' \nu''' p'''}
    \, _{ij}< K' \nu' p' | [V_{ij}-V_{hc}(\rho)] | K'' \nu'' p'' >_{ij}
    \times	\nonumber \\ & & \hskip 2truecm
    A_{K'' \nu'' p'', K''' \nu''' p'''}^{-1}
    \sum_{kl} \, _{ij}< K''' \nu''' p''' | K \nu p >_{kl}
    R_{K \nu p}(\rho) \; .
\end{eqnarray}
Using the hyperspherical coordinates, and the following definitions
\begin{eqnarray}
  T_{K'\nu'p',K\nu p} &=&
        _{ij}< K' \nu' p' | \hat K | K \nu p >_{ij}\,, \cr
  V^r_{K'\nu'p',K\nu p}(\rho) &=&
        _{ij}< K' \nu' p' | V_{ij}-V_{hc} | K \nu p >_{ij}\,,\cr
  W^{[3]}_{K'\nu'p',K\nu p} &=& \sum_{K'' \nu'' p''}
        A_{K' \nu' p', K'' \nu'' p''}^{-1} \,
        _{ij}< K'' \nu'' p'' | K \nu p >_{ik}\,,
        \hskip 2truecm k \neq i,j \,,\cr
  W^{[4]}_{K'\nu'p',K\nu p} &=& \sum_{K'' \nu'' p''}
        A_{K' \nu' p', K'' \nu'' p''}^{-1} \,
        _{ij}< K'' \nu'' p'' | K \nu p >_{kl}\,,
        \hskip 2truecm k,l \neq i,j \,,
\end{eqnarray}
Eq. (\ref{Fad-A}) can be expressed as an hyperradial equation
\begin{eqnarray} \label{radial}
  \sum_{\mu' } \Big[ &-& \frac{1}{2} A_{\mu \mu'}
  \big( \frac{d^2}{d\rho^2} + \frac{3 A -4}{\rho}\frac{d}{d\rho}
  \big) +\frac{1}{2}\frac{T_{\mu \mu'}}{\rho^2}
   +\frac{A(A-1)}{2}V_{hc}(\rho)A_{\mu \mu'}
   -E \, A_{\mu \mu'} \Big] R_{\mu'}(\rho) \cr
  &=&-\sum_{\mu' \mu''} V^r_{\mu \mu'} (\rho)
  \left( \delta_{\mu'\, \mu''}
  +2(A-2)W^{[3]}_{\mu', \mu''}
  +\frac{(A-2)(A-3)}{2} W^{[4]}_{\mu',\mu''}
  \right)  R_{\mu''}(\rho) \,.
\end{eqnarray}
Here $\mu \equiv (K \nu p)$, and $\hat K$ is the generalized,
hyperspherical, angular momentum operator.
This hyperradial equation can be solved by expanding the
hyperradial function $R_{\mu}(\rho)$ into basis functions
$\phi_{n_{\rho}}(\rho)$
\begin{equation}
   R_{\mu}(\rho) = \sum_{n_{\rho}}
                   C_{n_{\rho}}^{\mu} \phi_{n_{\rho}}(\rho)
\end{equation}
containing the Laguerre polynomials $L_{n_{\rho}}^a$
\begin{equation} \label{Frho}
   \phi_{n_{\rho}}(\rho) = \sqrt{n_{\rho}! \over (n_{\rho}+a)!} \,\,
    b^{-3({\it A}-1)\over 2}
    \left(\frac{\rho}{b}\right)^{{\it a}-(3 {\it A}-4)\over 2}
    L_{n_{\rho}}^a\left({\rho\over b}\right)\,
    \exp\left[-{\rho\over 2b}\right]\,.
\end{equation}
Multiplying the right and left hand sides by $\phi^\star_{n_{\rho}}(\rho)$
and performing the hyperradial integration one remains with
the generalized eigenvalue problem
\begin{equation}
\sum_{m} C_m ({\bf H}_{m n}-E {\bf M}_{m n})=0\,,
\end{equation}
where $m$ and $n$ stand for the sets ${K,\nu,p,n_{\rho}}$ and
${K',\nu',p',n_{\rho}'}$, respectively,
and the matrix elements ${\bf H}_{m n}$ and ${\bf M}_{m n}$ are
given by the corresponding hyperradial integrals.

Eq. (\ref{radial}) represents the main equation of the present method.
The following remarks have to be stressed in order here:

\begin{enumerate}

\item the 3(A-1)--dimensional Schr\"odinger equation has been reduced to
      a $3\times 4$ dimensional integro--differential equation, due to the
      presence of at most four--body terms ($W^{[4]}_{K'\nu'K\nu}$);

\item to the extent the convergence in the expansions is reached the result
      for $\Psi=\sum_{i<j} \Psi_{ij}$ can be considered
      as an "exact" ground state solution of the many--body
      Schr\"odinger equation.
\end{enumerate}
\section{CONSTRUCTION OF THE BASIS FUNCTIONS}\label{basis}

Technical difficulties in calculating the matrix elements
${\bf H}_{m n}$ and ${\bf M}_{m n}$ are encountered in the construction
of the  basis functions $H_{K\,\nu} (\Omega,s_i,t_i)$ for increasing numbers of
particles and hyperspherical angular momentum $K$, as well as in the
calculation of the two--, three-- and four--body terms that appear in
${\bf H}_{m n}$ and ${\bf M}_{m n}$.
The basis functions $H_{K\,\nu} (\Omega,s_i,t_i)$ must be totally
antisymmetric functions in total space (coordinate, spin, isospin).
Therefore, one needs efficient
algorithms to construct convenient expressions for the symmetrized
spin--isospin basis functions as well as symmetrized hyperspherical
wave functions.

\subsection{Spin--isospin states with arbitrary permutational symmetry}

The construction of A--particle symmetrized spin--isospin
functions as well as hyperspherical functions
can be done recursively \cite{Akiva88}.
Each  A--particle spin--isospin function of a well defined
irrep of the symmetric group
$\Gamma_A$ and spin--isospin quantum numbers $S_A, T_A$,
is written as a linear combination of spin--isospin coupled products
of an (A-1)--particle wave function of well defined irrep $\Gamma_{A-1}$
and spin--isospin quantum numbers $S_{A-1}, T_{A-1} $,
and the A'th particle wave function with spin--isospin quantum numbers
$s_A, t_A$.
The coefficients of this
linear combination are a sort of coefficients of fractional parentage.
The subspace of functions with good quantum numbers $S_A, T_A$ and
$\Gamma_{A-1}$ is an invariant subspace for the transposition class--sum
operator of the symmetric group ${\cal S}_A$. The eigenvalues that are obtained
after the diagonalization of the class--sum operator identify
the irreps of the symmetric group uniquely. The eigenvectors are the
spin--isospin coefficients of fractional parentage (stcfps)
\begin{eqnarray} \label{stcfps}
& & {\cal X}_{S_A {S^z}_A T_A {T^z}_A \,
    {\Gamma}_{A} Y_{A-1}  \, \alpha^{ST}_A} (s_i t_i)
 =  \sum_{S_{A-1} T_{A-1} \alpha^{TS}_{A-1} }
\nonumber \\ & & \hskip 2truecm
   [ (S_{A-1}; s_A)S_A (T_{A-1} ; t_A)T_A \Gamma_{A-1}\alpha^{ST}_{A-1}
   |\} S_A T_A {\Gamma}_{A} \alpha^{ST}_A ] \times
\nonumber \\ & & \hskip 2truecm
  \big[ {\cal X}_{S_{A-1} {S^z}_{A-1} T_{A-1} {T^z}_{A-1} \,
  {\Gamma}_{A-1} Y_{A-2}
  \, \alpha^{ST}_{A-1} } \otimes s_{A} t_A \big]^{S_A {S^z}_A T_A {T^z}_A}\,.
\end{eqnarray}
The sum over $S_{A-1}$ and $T_{A-1}$ in Eq. (\ref{stcfps}) is subject to the
usual angular momentum coupling rules.

\subsection{Symmetrized hyperspherical functions}

In an analogous way, one can use the recursive methods
developed in the last few years  \cite{{Akiva94},{Nir98b}} for
constructing hyperspherical functions that belong to well defined
irreps of the symmetric group.
In this methods the reduction problem
${\cal O}_{3A-3} \supset {\cal O}_3 \otimes {\cal S}_A$ is solved and one obtains
hyperspherical functions
which belong to irreps of the symmetric group ${\cal S}_A$ and have
good angular momentum and hyperspherical angular momentum quantum numbers.
These functions are expressed in terms of ${\cal S}_{A-1}$ hyperspherical states,
that are coupled, via the ``tree" method
 ~\cite{Vilenkin66},
with the appropriate single--Jacobi coordinate hyperspherical functions,
into coupled $A$--particle states with the desired
angular momentum and hyperspherical angular momentum quantum numbers.
Actually, these coupled states yield an invariant subspace with
respect to ${\cal S}_A$. The transposition class--sum of the symmetric
group (the second Casimir operator) is diagonalized within this subspace.
The eigenvalues that are obtained after the diagonalization
identify the irreps of the symmetric group uniquely, and
the eigenvectors are the hyperspherical coefficients of fractional
parentage (hscfps).

Constructing basis functions in such a recursive way makes the evaluation
of any two body operator easy, since only the matrix element of the two--body
operator between the last two particles needs to be calculated.
A further improvement in the efficiency of the algorithm
can be reached if one uses reversed order Jacobi coordinates for
the construction of the hyperspherical functions.
Normally Jacobi coordinates are defined so that the first
one is the relative distance between particles 1 and 2, while the last
is the distance between the A'th particle and the center of mass of the
A-1 particle system. This implies that in order to calculate the matrix
element of a two--body operator, depending typically
on the interparticle distance between the A'th and (A-1)'th particles (see
Section \ref{matrix}),
one needs to rotate the last two Jacobi coordinates by a proper angle so
that one of them represents the interparticle distance. This can be done
by using the Raynal Revai, the T-coefficients and the 6-j coefficients.
Constructing the Jacobi coordinates in reverse order, would
simplify the calculation of two-- and three--body matrix elements as no further
rotation is needed \cite{Nir98b}.

Recently an alternative way to construct the symmetrized
hyperspherical functions has been
proposed \cite{Nir97}. This method consists in introducing as an
intermediate subgroup the orthogonal group of kinematic rotations
${\cal O}_{A-1}$,
i.e. one uses the group chain
${\cal O}_{3A-3}\supset {\cal O}_3 \otimes {\cal O}_{A-1}\supset
{\cal O}_3 \otimes {\cal S}_A$.
The introduction of the kinematic group turns out to be of particular
importance for increasing number of particles and for large values of $K$
as it results in a significant reduction in the number of hscfps and in
the computation time.
Another benefit in using this method is the realization of kinematic
rotations through the representation matrices of the group ${\cal O}_{A-1}$
thus avoiding the use of the Raynal--Revai and the T--coefficients
in calculating matrix elements of two--body operators depending on the
interparticle distance.
In the present calculation we have adopted this procedure, even if it is
less efficient than defining Jacobi coordinates in reversed order. The
reason is that in this way one can easily extend the present method
to solve the Schr\"odinger--like equation with source, necessary for
the calculation of the Lorentz integral transforms of response functions ~\cite{ELO94}.
The presence of the source annuls the advantages of the alternative set of
Jacobi coordinates.

Summarizing the procedure for the construction of the symmetrized
hyperspherical harmonics one can split it into two steps.
In the first step one constructs
hyperspherical functions with good orthogonal symmetry, i.e. basis
functions with good quantum numbers $K_N, L_N, M_N$ that belong
to a well defined Gel'fand--Zetlin pattern
$\mbox{\boldmath $\Lambda$}_{N}=
	 (\mbox{\boldmath $\lambda$}_{N},
          \mbox{\boldmath $\lambda$}_{N-1},\ldots,
          \mbox{\boldmath $\lambda$}_{2})$.
The orthogonal group irreps $\mbox {\boldmath $\lambda$}_j$,($j=2,3,\ldots $)
are characterized by the integer or half integer
numbers
$\lambda_{j,1}$, $\lambda_{j,2}$, $ \ldots $ , $\lambda_{j,k}$,
where $k=[\frac{j}{2}]$.
In Ref. \cite{Nir97} it was pointed out that for hyperspherical
functions these
numbers are always integers and, in addition, there are at
most three non--zero
values in the irrep \mbox{\boldmath $\lambda$}$_j$ for $j \geq 6 $ i.e.,
$\lambda_{j,1}$, $\lambda_{j,2}$, and $\lambda_{j,3}$.
The second step is to reduce each irrep of the orthogonal group
${\cal O}_N$ into
irreps of the symmetry group ${\cal S}_{N+1}$.
These two steps are carried out using the recursive method developed
in Ref. \cite{Nir97}.
At the end, the symmetrized A--particle,
$N$--Jacobi coordinate, hyperspherical states can be expressed in terms
of $(A-1)$--particle states coupled to the $A$'th particle state by using
two new types of coefficients of
fractional parentage, namely the orthogonal--hyperspherical cfps (ohscfps)
for the construction of hyperspherical functions with good orthogonal
symmetry and the orthogonal cfps (ocfps) for the reduction
${\cal O}_{A-1} \downarrow {\cal S}_{A}$.
One has
\begin{eqnarray} \label{symhh}
\label{twocfp} & &
{\cal Y}_{
  K_{N} L_{N} M_{N} \mbox {\boldmath $\lambda$}_{N} \beta^K_{N}
  \Gamma_A Y_{A-1}  \beta^{\lambda}_A} (\Omega)
\nonumber \\ & &
= \sum_{ {\tiny \mbox {\boldmath $\lambda$}_{N-1}} \beta^{\lambda}_{A-1} }
 [ \left( \mbox {\boldmath $\lambda$}_{N-1} \Gamma_{A-1}
 \beta^{\lambda}_{A-1} \right)
   \mbox{\boldmath $\lambda$}_{N}
  |\} \mbox {\boldmath $\lambda$}_{N} \Gamma_A \beta^{\lambda}_A ]\times
\nonumber \\ & &
\sum_{K_{N-1} L_{N-1} \beta^K_{N-1}  \ell_{N} }
  \left[  \big( K_{N-1} L_{N-1} \mbox {\boldmath $\lambda$}_{N-1}
  \beta^K_{N-1};\ell_{N} \big)
  K_{N} L_{N} |\} K_{N} L_{N} \mbox {\boldmath $\lambda$}_{N}
\beta^K_{N} \right]
 \nonumber \\ & & \; \; \;
{\cal Y}_{
   \big( K_{N-1} L_{N-1}  \mbox {\boldmath $\lambda$}_{N-1} \beta^K_{N-1}
   \Gamma_{A-1} Y_{A-2}   \beta^{\lambda}_{A-1} ; \ell_{N} \big)
   K_{N} L_{N} M_{N} }(\Omega) \; \; .
\end{eqnarray}
Note that the order of the summation
in this equation is important. We have to start from the sum over the
ocfps
\begin{equation}
[ \left( \mbox {\boldmath $\lambda$}_{N-1} \Gamma_{A-1}
 \beta^{\lambda}_{A-1} \right)  \mbox{\boldmath $\lambda$}_{N}
  |\} \mbox {\boldmath $\lambda$}_{N} \Gamma_A \beta^{\lambda}_A ] \, ,
\nonumber
\end{equation}
to determine the irrep $\mbox{\boldmath $\lambda$}_{N-1}$ of ${\cal O}_{N-1}$
and then we can sum over the ohscfps
\begin{equation} \label{cfpON}
\left[  \big( K_{N-1} L_{N-1} \mbox {\boldmath $\lambda$}_{N-1}
  \beta^K_{N-1};\ell_{N} \big)  K_{N} L_{N} |\} K_{N} L_{N}
  \mbox {\boldmath $\lambda$}_{N} \beta^K_{N} \right] \, ,
\nonumber
\end{equation}
where $(\beta^K_N$ and $\beta^{\lambda}_A)$ are the degeneracy removing
labels ($(\beta^K_N,\beta^{\lambda}_A)\equiv\alpha^K_N$).

\section{Calculation of the matrix elements}\label{matrix}

The matrix elements of any two--body operator
\begin{equation} \label{O_ij}
   O_{ij} = O^R(\vec{r}_{ij}) O^{ST}(s_i,s_j,t_i,t_j)
\end{equation}
between the fermionic, antisymmetric, hyperspherical basis functions of
Eq. (\ref{H}), can be written as a sum of a spatial term multiplied by
a spin--isospin term,
\begin{eqnarray} \label{2bme}
 & & < K \nu | O_{ij} | K' \nu' > =
      \sum_{\Gamma_{A-1}\Gamma'_{A-1}\Gamma_{A-2}}
      \prod (\mbox{ ISF })
      \times \cr & & \hskip 2truecm
     < S_A {S^z}_A T_A {T^z}_A \, \tilde{\Gamma}_{A} \tilde{Y}_{A-1}
        \, \alpha^{ST}_A |   O^{ST}(s_i,s_j,t_i,t_j)
      | S'_A {S'^z}_A T'_A {T'^z}_A \, \tilde{\Gamma}'_{A} \tilde{Y}'_{A-1}
        \, \alpha'^{ST}_A >
      \times \cr & & \hskip 0.5truecm
      < K_{N} L_{N} M_{N} \mbox {\boldmath $\lambda$}_{N}
       \beta^K_{N} \Gamma_A Y_{A-1}  \beta^{\lambda}_A |
        O^R(\vec{r}_{ij})
      | K'_{N} L'_{N} M'_{N} \mbox {\boldmath $\lambda$}'_{N}
       \beta'^K_{N} \Gamma'_A Y'_{A-1}  \beta'^{\lambda}_A > \,,
\end{eqnarray}
where the sum runs over all the ${\cal S}_{A-2}, {\cal S}_{A-1}$ irreps
$\Gamma_{A-1}\Gamma'_{A-1}\Gamma_{A-2}$ such that
$\Gamma_{A-2}\in\Gamma_{A-1}\in\Gamma_{A}$ and
$\Gamma_{A-2}\in\Gamma'_{A-1}\in\Gamma'_{A}$.
The factor preceding the spatial and spin--isospin matrix element is
a product of the inner product symmetric group isoscalar factors for
the antisymmetric representation and is given by
\begin{equation}
 \prod( \mbox{ ISF }) =
    \Lambda_{\Gamma_{A}\Gamma_{A-1}}
         \sqrt{\frac{|\Gamma_{A-1}|}{|\Gamma_{A}|}}
    \Lambda_{\Gamma_{A-1}\Gamma_{A-2}}
         \sqrt{\frac{|\Gamma_{A-2}|}{|\Gamma_{A-1}|}}
    \Lambda_{\Gamma'_{A}\Gamma'_{A-1}}
         \sqrt{\frac{|\Gamma'_{A-1}|}{|\Gamma'_{A}|}}
    \Lambda_{\Gamma'_{A-1}\Gamma'_{A-2}}
         \sqrt{\frac{|\Gamma'_{A-2}|}{|\Gamma'_{A-1}|}}\, .
\end{equation}
The phase factor
 $\Lambda_{\Gamma_{A}  \Gamma_{A-1}}$ \cite{Akiva88b}
is positive (negative) when the number of boxes in $\Gamma_{A}$
below the row of the $A$'th particle is even (odd).

Using the permutation symmetry of the HH states the matrix elements of
any two--body operator, $O_{ij}$, are equal to the matrix elements of
$O_{A \,,A-1}$. The spatial matrix elements are then calculated
using the last two generations of cfps, referring to the
constructions $ A-2 \longrightarrow A-1$ and
$ A-1 \longrightarrow A $,
and the appropriate kinematical rotations which
reduce the calculation of any
two--body operator to sum over one dimensional integrals.
The formal derivation of the calculation
of two--body operator matrix elements between
the symmetrized hyperspherical harmonics of Eq. (\ref{symhh}) has been
presented in rather a detailed manner by Barnea and Novoselsky
\cite{Nir98a}. Summarizing their results one has
\begin{eqnarray} \label{DVN}
 & & < K_{N} L_{N} M_{N} \mbox {\boldmath $\lambda$}_{N}
       \beta^K_{N} \Gamma_A Y_{A-1}  \beta^{\lambda}_A |
        O^R(\vec{r}_{A,A-1})
      | K'_{N} L'_{N} M'_{N} \mbox {\boldmath $\lambda$}'_{N}
       \beta'^K_{N} \Gamma'_A Y'_{A-1}  \beta'^{\lambda}_A > =
\cr & & \hskip 2truecm
	\delta_{\Gamma_{A-2}, \Gamma'_{A-2}}
	\sum_{\mbox {\boldmath $\lambda$}_{N-1} \beta^{\lambda}_{A-1}}
	\sum_{\mbox {\boldmath $\lambda$}'_{N-1}
              \beta^{\lambda \, \prime}_{A-1}}
	\sum_{\mbox {\boldmath $\lambda$}_{N-2} \beta^{\lambda}_{A-2}}
\cr & & \hskip 2truecm
    [ \left( \mbox {\boldmath $\lambda$}_{N-1} \Gamma_{A-1}
      \beta^{\lambda}_{A-1} \right)  \mbox{\boldmath $\lambda$}_{N}
      |\} \mbox {\boldmath $\lambda$}_{N} \Gamma_A \beta^{\lambda}_A ]^*
    [ \left( \mbox {\boldmath $\lambda$}'_{N-1} \Gamma'_{A-1}
      \beta^{\lambda \, \prime}_{A-1}\right) \mbox{\boldmath $\lambda$}'_{N}
    |\} \mbox {\boldmath $\lambda$}'_{N} \Gamma'_A \beta^{\lambda\,\prime}_A ]
\cr & & \hskip 1truecm \times
    [ \left( \mbox {\boldmath $\lambda$}_{N-2} \Gamma_{A-2}
      \beta^{\lambda}_{A-2} \right)  \mbox{\boldmath $\lambda$}_{N-1}
    |\} \mbox{\boldmath $\lambda$}_{N-1}\Gamma_{A-1}\beta^{\lambda}_{A-1}]^*
    [ \left( \mbox {\boldmath $\lambda$}_{N-2} \Gamma_{A-2}
      \beta^{\lambda}_{A-2} \right)  \mbox{\boldmath $\lambda$}'_{N-1}
    |\} \mbox{\boldmath $\lambda$}'_{N-1}\Gamma'_{A-1}
        \beta^{\lambda\,\prime}_{A-1}]
\cr & & \hskip 2truecm \times
    \sum_{\mbox {\boldmath $\lambda$}''_{N-1}}
    \mbox { $D$ }_{ \tilde{\bf \Lambda}_{N-1}^{ \prime \prime }
                    \tilde{\bf \Lambda}_{N-1} }
                 ^{  \mbox {\boldmath $\lambda$}_{N} \; * }
                  (g_{A,A-1})
    \mbox { $D$ }_{ \tilde{\bf \Lambda}_{N-1}^{ \prime \prime }
                    \tilde{\bf \Lambda}'_{N-1} }
                 ^{  \mbox {\boldmath $\lambda$}'_{N} }
                  (g_{A,A-1})
\cr & & \hskip 1truecm \times
    < K_{N} L_{N} M_{N} \mbox {\boldmath $\lambda$}_{N} \beta^K_{N}
    \tilde{\bf \Lambda}^{\prime \prime} _{N-1}
    Y_{A-2}  \beta^{\lambda}_{A-2})
    |  O^R (\sqrt{2} \vec{\xi}_{N} )
    | K'_{N} L'_{N} M'_{N} \mbox {\boldmath $\lambda$}'_{N}
    \beta^{K \, \prime}_{N} \tilde{\bf \Lambda}^{ \prime \prime}_{N-1}
    Y_{A-2}  \beta^{\lambda}_{A-2}) > \; .
\end{eqnarray}
Here, for simplicity, we use the symbol ${\bf \Lambda}_{N-1} $
to denote the irreps
$ [\mbox {\boldmath $\lambda$}_{N-1}, \mbox {\boldmath $\lambda$}_{N-2}]$.
Note that the sum in Eq. (\ref{DVN}) over the irreps
$ \mbox {\boldmath $\lambda$}^{\prime \prime}_{N-1} $ is restricted to those
states which belong both to the irrep $ \mbox {\boldmath $\lambda$}_{N}$
and $ \mbox {\boldmath $\lambda$}'_{N}$.
The kinematical ``relative coordinate'' rotation applied to the Jacobi coordinates
\begin{eqnarray} \label{rotation1}
\vec{\xi}_{N-1}  & = & \sqrt{\frac{2 (A-2)}{A}}
                       \big( \frac{\vec{r}_{A-1}+\vec{r}_{A}}{2}
			  -\frac{1}{A-2}\sum_{i=1}^{A-2} \vec{r}_i \big)
\nonumber \\
\vec{\xi}_{N}    & = & \sqrt{\frac{1}{2}}
                       \big( \vec{r}_{A}-\vec{r}_{A-1} \big)
                 \; ,
\end{eqnarray}
and given by
\begin{eqnarray} \label{rotation2}
\vec{\xi}_{N-1}  & = & \sqrt{\frac{A}{2(A-1) } }\vec{\eta}_{N-1}  +
                       \sqrt{\frac{A-2}{2(A-1)}}\vec{\eta}_{N}
\nonumber \\
\vec{\xi}_{N}    & = & -\sqrt{\frac{A-2}{2(A-1)}}\vec{\eta}_{N-1}
                       +\sqrt{\frac{A}{2(A-1) } }\vec{\eta}_{N}
                 \; ,
\end{eqnarray}
is realized in Eq. (\ref{DVN}) by the generalized, ${\cal O }_N$, Wigner
$D$ functions, where $g_{A,A-1}$ is the group element that corresponds to
the rotation (\ref{rotation2}), i.e. a rotation by an angle
$\gamma=\arcsin(\sqrt{\frac{A-2}{2(A-1)}})$ in the
$\vec{\eta}_{N},\vec{\eta}_{N-1}$ plane.

Since the two--body operator in Eq. (\ref{DVN}) depends only on the
coordinate $ \vec{\xi}_{N}$ we should
separate the hyperspherical functions related to this coordinate
in the bra and the ket states. So, the last step in the derivation
consists of using the ohscfps introduced in
Eq. (\ref{cfpON}). Then, the last term on the rhs of Eq.  (\ref{DVN}) is
\begin{eqnarray} \label{lasteq}
&  &     < K_{N} L_{N} M_{N} \mbox {\boldmath $\lambda$}_{N} \beta^K_{N}
    \tilde{\bf \Lambda}^{\prime \prime} _{N-1}
    Y_{A-2}  \beta^{\lambda}_{A-2})
    |  O^R (\sqrt{2} \vec{\xi}_{N} )
    | K'_{N} L'_{N} M'_{N} \mbox {\boldmath $\lambda$}'_{N}
    \beta^{K \, \prime}_{N} \tilde{\bf \Lambda}^{ \prime \prime}_{N-1}
    Y_{A-2}  \beta^{\lambda\,\prime}_{A-2}) > =
\nonumber \\ & &
   \sum_{K_{N-1} L_{N-1} \beta^K_{N-1}  \ell_{N}, \ell'_{N} }
  \left[  \big( K_{N-1} L_{N-1} \mbox {\boldmath $\lambda$}_{N-1}
  \beta^K_{N-1};\ell_{N} \big)  K_{N} L_{N} | \}
  K_{N} L_{N} \mbox {\boldmath $\lambda$}_{N} \beta^K_{N} \right]
  \nonumber \\ & & \hskip 2truecm
\times \left[  \big( K_{N-1} L_{N-1} \mbox {\boldmath $\lambda$}_{N-1}
  \beta^K_{N-1};\ell'_{N} \big)  K'_{N} L'_{N} | \}
  K'_{N} L'_{N} \mbox {\boldmath $\lambda$}'_{N} \beta^{K \, \prime}_{N}
  \right]
  \nonumber \\ & & \hskip 1truecm \times
  < \big( K_{N-1} L_{N-1} \ldots; \ell_{N} \big) K_{N} L_{N} M_{N}  |
         O^R(\sqrt{2}\vec{\xi}_{N})
  | \big( K_{N-1} L_{N-1} \ldots  ; \ell'_{N} \big)
    K'_{N} L'_{N} M'_{N}  > \; \; .
\end{eqnarray}
This matrix element can be calculated for any given two--body operator.
For scalar operator $O^R(r_{A,A-1})$, one obtains
\begin{eqnarray} \label{scalar2bme}
\lefteqn{ O^R_{K_{N} K'_{N}; K_{N-1} \ell_N}(\rho) \equiv }
    \nonumber \\ & &
  < \big( K_{N-1} L_{N-1} \ldots; \ell_{N} \big) K_{N} L_{N} M_{N}  |
         O^R(\sqrt{2}{\xi}_{N})
  | \big( K_{N-1} L_{N-1} \ldots  ; \ell'_{N} \big)
    K'_{N} L'_{N} M'_{N}  > =
    \cr & & \hskip 1truecm
    \delta_{\ell_N, \ell'_N}
    {\cal N}_{N}( K_{N} ; \ell_{N} K_{N-1} )
    {\cal N}_{N}( K'_{N} ; \ell_{N} K_{N-1} )
  \; \int _0^{\frac{\pi}{2}} d \theta
  \sin^{2 \ell_{N}+2}(\theta) \cos^{2 K_{N-1}+3N-4}(\theta)
\cr & & \hskip 2truecm \times
  P_{n}
   ^{(\ell_{N}+\frac{1}{2},K_{N-1}+\frac{3N-5}{2}) }(\cos{2\theta} )
  P_{n'}
   ^{(\ell_{N}+\frac{1}{2},K_{N-1}+\frac{3N-5}{2}) }(\cos{2\theta} )
   O^R(\sqrt{2} \rho \sin \theta)  \, ,
\end{eqnarray}
with $n=\frac{K_N-\ell_N-K_{N-1}}{2}$,
$n'=\frac{K'_N-\ell_N-K_{N-1}}{2}$ and where
\begin{equation} \label{norN}
  \bbox {\cal N}_N(K_N ; \ell_N K_{N-1}) = \left[ \frac
  {(2K_N+3N-2)n!\Gamma(n+K_{N-1}+\ell_N+\frac{3N-2}{2})}
  {\Gamma(n+\ell_N+\frac{3}{2})
  \Gamma(n+K_{N-1}+\frac{3N-3}{2})}
  \right]^{\frac{1}{2} }
\end{equation}
is a normalization constant, and $P_n^{(a,b)}$ are the Jacobi polynomials.

The spin--isospin part of the two--body operator is calculated in a similar
way, first we use the last two generations of the spin--isospin cfps (
Eq. (\ref{stcfps})) to get the explicit dependence of the state
$| S_A {S^z}_A T_A {T^z}_A \, \tilde{\Gamma}_{A} \tilde{Y}_{A-1}
        \, \alpha^{ST}_A >$
on the spin and the isospin of the last two--particles,
\begin{eqnarray} \label{stme1}
& &< S_A {S^z}_A T_A {T^z}_A \,{\Gamma}_{A} {Y}_{A-1} \, \alpha^{ST}_A |
    O^{ST}(s_A,s_{A-1},t_A,t_{A-1})
   | S'_A {S'^z}_A T'_A {T'^z}_A \,{\Gamma}'_{A} {Y}'_{A-1}\,\alpha'^{ST}_A >
   = \cr & & \hskip 1truecm
   \delta_{Y_{A-2}, Y'_{A-2}}
   \sum_{S_{A-1}T_{A-1}\alpha^{ST}_{A-1}}
   \sum_{S'_{A-1}T'_{A-1}\alpha'^{ST}_{A-1}}
   [ (S_{A-1}; s_A)S_A (T_{A-1} ; t_A)T_A \Gamma_{A-1}\alpha^{ST}_{A-1}
   |\} S_A T_A {\Gamma}_{A} \alpha^{ST}_A ]
   \cr & & \hskip 1truecm \times
   [(S'_{A-1}; s_A)S'_A (T'_{A-1} ; t_A)T'_A \Gamma'_{A-1}\alpha'^{ST}_{A-1}
   |\} S'_A T'_A {\Gamma}'_{A} \alpha'^{ST}_A ]
   \cr & & \hskip 1truecm \times
   \sum_{S_{A-2}T_{A-2}\alpha^{ST}_{A-2}}
   [ (S_{A-2}; s_A)S_{A-1} (T_{A-2} ; t_A)T_{A-1}\Gamma_{A-2}
   \alpha^{ST}_{A-2} |\} S_{A-1} T_{A-1} {\Gamma}_{A-1} \alpha^{ST}_{A-1} ]
   \cr & & \hskip 1truecm \times
   [ (S_{A-2}; s_A)S'_{A-1} (T_{A-2} ; t_A)T'_{A-1}\Gamma_{A-2}
   \alpha^{ST}_{A-2} |\} S'_{A-1} T'_{A-1} {\Gamma}'_{A-1} \alpha'^{ST}_{A-1} ]
   \\ & & \hskip 1truecm \times
   < S_A T_A S_{A-1} T_{A-1} S_{A-2} T_{A-2} \ldots |
   O^{ST}(s_A,s_{A-1},t_A,t_{A-1})
   | S'_A T'_A S'_{A-1} T'_{A-1} S'_{A-2} T'_{A-2} \ldots > \nonumber
   \,.
\end{eqnarray}
Then the matrix element on the rhs of (\ref{stme1}) is calculated
using the $6j$ symbols.

Special attention should be payed to the matrix elements of the generalized
angular momentum operator $\hat{K}$ as it contains the hyperangular part
of the Laplace operator. Using the definition in Eq. (\ref{defME}) we
see that
\begin{equation}
  T_{K'\nu'p',K\nu p} =
        < K' \nu' | \chi_{p'}(z_{ij})^* \hat K \chi_p(z_{ij})
    | K \nu > \, .
\end{equation}
This matrix element can be easily calculated once we know the action
of $\hat K$ on $\chi_p(z_{ij})$. Since $\hat K$ is invariant under
kinematical rotations we can always write it using a set of Jacobi coordinates
such that $\vec\xi_N=\frac{1}{\sqrt{2}}(\vec{ r}_i-\vec{ r}_j)$, which leads
to the following expression
\begin{eqnarray}
& &< K' \nu' | \chi_{p'}(z_{ij})^* \hat K \chi_p(z_{ij}) | K \nu > =
   \cr & & \hskip 2truecm
   < K' \nu' | \chi_{p'}(z_{ij})^* K(K+3N-2) \chi_p(z_{ij})
    - \chi_{p'}(z_{ij})^* (1-z^2_{ij})\frac{d^2 \chi_p(z_{ij})}{d z_{ij}^2}
   \cr & & \hskip 3truecm
    + \chi_{p'}(z_{ij})^* \{ z_{ij} + \frac{-2+(3N-2)z^2_{ij}}{z_{ij}} \}
      \frac{d \chi_p(z_{ij})}{d z_{ij}}
   \cr & & \hskip 3truecm
    - 2 \chi_{p'}(z_{ij})^*
     (1-z^2_{ij})\frac{d \chi_p(z_{ij})}{d z_{ij}} \frac{d}{d z_{ij}}
     | K \nu > \,.
\end{eqnarray}
This is just the matrix element of a sum of two body operators which can
be calculated as explained above.

A bit more problematic are the three-- and four--body integrals
$W^{[3]}_{K'\nu'\,p'\, K \nu\,p}$ and
$W^{[4]}_{K'\nu'\,p' K\nu\,p}$. In order to calculate the three-- (four--)
body term we must first use three (four) generations of the ocfps in order
to express the matrix elements in terms of the orthogonal symmetry
adapted hyperspherical functions. Then we can use the proper rotations
and reduce the integrals to six--dimensional integrals that depend
only on the last two Jacobi coordinates. The explicit dependence of
the HH functions on the Jacobi coordinates is then revealed using the
last two generations of the ohscfps, Eq. (\ref{cfpON}).
For $L=0$ states, these
integrals can further be reduced to three--dimensional integrations.

\section{DISCUSSION OF THE RESULTS}\label{results}

The present method becomes more and more complex as the number
of fermions increases. Therefore in this work we consider
only central $NN$ potentials. We present results
for 6, 8 and 12 nucleons interacting via the Volkov ~\cite{Volkov} (VV),
the Afnan--Tang ~\cite{Afnan} (S3), the modified S3 potential
~\cite{ModAfnan} (MS3), the Brink--Boeker ~\cite{Brink} (B1)
and the Malfliet--Tjon ~\cite{Malfliet} potentials (MT--I/III and MT--V).
The accuracy and convergence of the method have been investigated
for 6 and 8 nucleons.

The first step in our numerical study was to determine the
hyperradial working point, i.e. the parameters $a$ and $b$ of the Laguerre
polynomials in Eq. (\ref{Frho}), the number of radial grid points and
their location.
After some preliminary tests we decided to use 15 hyperradial grid
points and 15 hyperradial functions
with the parameters $b=0.25$ fm and $a=12$. The grid points
where chosen as the abscissas for a Gauss--Laguerre integration.

The next step in setting our working point was to check the convergence
of the two--body correlation terms. As an example we studied the $^8$Be
system with the hard core MT--V potential, using a single
HH function, with $K=4$ and expanding the correlation function into
Jacobi polynomials.
As can be seen from Fig. \ref{fig:corr} the binding energy converges
very fast with increasing number $n_x$ of polynomials .
In fact, the energy difference between
the calculations with 7 and 8 correlations terms is about 0.02 MeV.
As a consequence we used 8 terms in the expansion of the correlation
function.

The effect of introducing the hypercentral potential $V_{hc}(\rho)$ is
presented in Table \ref{tb:Vhc} for $^6$Li with the VV
 potential. As can be seen the hypercentral potential accelerates
the convergence of the HH expansion. With the lowest $K_{max}$ of 2
one already obtains a larger binding energy than with $K_{max}=6$
without hypercentral potential. Thus the introduction of $V_{hc}$
is certainly advantageous.

Our numerical results with the parameters described above are presented
in Table \ref{tb:n6} for the $6$ nucleon system and in Table \ref{tb:n8}
for the $8$ nucleon system.

For $^6$Li the calculations include only one line
irreps of the kinematical group ${\cal O}_5$ and
the irreps [42] and [33] of the permutation group.
We compare our values with recent accurate variational
results  ~\cite{Varga95} available for some of these central potentials.
One can notice that with the VV potential one reaches convergence
faster than with other potentials. The result for the binding energy
starts oscillating around the asymptotic value. It also compares nicely
to the variational result of Ref. ~\cite{Varga95}. The other potentials
show a tendency to convergence even if $K_{max}$ is not large enough to
reach it. The MT potentials seem to lead to a more rapid convergence
than S3 and B1. The differences may be due
both to the fact that the convergent value has not yet been reached and
to the missing irreps of the permutation group.
The two line irreps of the orthogonal group are of little importance.
This has been checked for
the 8 particle case, where they give rather small contributions.

The calculations for $^8$Be
include only the irrep (400) of the kinematical group ${\cal O}_7$  and
irreps with at most 3 rows of the permutation group. Here the
comparison is made with the TICI results ~\cite{GMNBPW96}. In all cases our
results for the binding energy are somewhat larger. They
show characteristics similar to the six--body case. Again one sees
that the VV potential
result presents small oscillations around the convergent value and
that the rather hard core MT potentials  seem to give
values closer to convergence than B1 or S3.
From the comparison between the TICI results with and without state dependent
correlations one can infer that for the MT potentials state independent
correlations already give rather satisfying results, while
state dependent correlations are more effective for
B1 and S3. Since our correlations are state independent one could expect such
a different convergence behavior as found in Table  \ref{tb:n8}

The results for $^{12}$C are presented in Table \ref{tb:n12}.
Here the calculations include a single HH state, $K_{max}=8$,
with the kinematical group ${\cal O}_{11}$
irrep $(4,4,0)$ and the ${\cal S}_{12}$ irrep $[444]$.
The comparison is made with TICI ~\cite{GMNBPW96}, IDEA ~\cite{Brizzi96},
and the variational result of Ref. ~\cite{BBG87}.
Again from the comparison with the TICI results one can see that S3 and B1
results are farer from convergence than the MT values and that
state dependent correlations play a similar role as discussed
for the $^8$Be case.

The IDEA results for 12 particles deviate somewhat from ours, although
they are all obtained using a single HH state. One possible explanation of
the difference can be attributed to the choice of the HH state which
is not unique.
On the other hand we have obtained excellent agreement with the IDEA results
for 16 bosons.

\section{CONCLUSIONS}\label{concl}

In this work we have formulated a general method
to calculate the wave functions of light systems up to considerably
large number of particle. This method combines the main ideas of
the HH expansion, the pair correlation ansatz and the Faddeev approach.
The actual application of it is made possible by the use of a
very efficient recursive algorithm to construct the antisymmetric A-particle
state containing hyperspherical harmonics. We have applied the method to
calculate the binding energies of 6, 8 and 12 nucleon systems with
central local potentials.
The results we have obtained are very encouraging. For some potentials
(VV, 6 particles and MTV, 8 particles) we have reached the
convergence region with $K_{max}=8$ which is the maximum value allowed
by our present computer facilities (workstations).
The $^6$Li result for the VV potential is slightly higher than the SVM result.
Even if in other cases we  have not yet reached the convergence in the HH expansion our results are close to the TICI results. For the eight--body case
they are higher for all the potentials where results were
available for a comparison.

The method  presented here for the solution of the few--body Schr\"odinger
equation can be easily extended to solve the Schr\"odinger--like
equation with a source, necessary for the application of the Lorentz integral
transform method. Work in this direction is in progress.
\begin{table}
\caption{\label{tb:Vhc} \em \\
         The effect of the hypercentral potential $V_{hc}$ on the
         binding energy of six--nucleon system $^6$Li,
         $(L,S)J^{\pi}=(0,1)1^{+}$, interacting via the VV
         potential. }
\begin{center}
\begin{tabular}{lcc}
  $K_{max}$ & \hspace{6mm} $V_{hc}=0$    \hspace{6mm}
	    & \hspace{6mm}
         $V_{hc}=\int d \Omega \, H_{K_{min} \nu}^{\dagger}
                 V_{1,2} H_{K_{min}\nu} $  \hspace{6mm} \\ \hline
 2 & 64.18 & 66.10  \\
 4 & 64.81 & 66.53  \\
 6 & 65.47 & 66.63  \\
\end{tabular}
\end{center}
\end{table}

\begin{table}
\caption{\label{tb:n6} \em \\
         Binding energies of six--nucleon system $^6$Li,
         $(L,S)J^{\pi}=(0,1)1^{+}$, interacting via various
         NN potentials. $N_{HH}$ represents the number of hyperspherical
         harmonic states.}
\begin{center}
\begin{tabular}{lcccccc}
  $K_{max}$ & $N_{HH}$ & \hspace{6mm} B1       \hspace{6mm}
		       & \hspace{6mm} MT--I/III \hspace{6mm}
		       & \hspace{6mm} MTV      \hspace{6mm}
		       & \hspace{6mm} S3       \hspace{6mm}
                       & \hspace{6mm} VV   \hspace{6mm} \\ \hline
 2 &  1  & 30.99 & 30.15 & 62.45 & 62.76 & 66.10  \\
 4 &  4  & 37.82 & 34.67 & 63.27 & 64.45 & 66.53  \\
 6 & 12  & 39.11 & 35.43 & 64.10 & 66.49 & 66.63  \\
 8 & 31  & 39.61 & 35.91 & 64.55 & 67.18 & 66.57  \\
\hline
SVM \cite{Varga95}  &
         &   -    &  -    & 66.30 & 70.65 & 66.25  \\
\end{tabular}
\end{center}
\end{table}

\begin{table}
\caption{\label{tb:n8} \em \\
         Binding energies of eight--nucleon system $^8$Be,
         $(L,S)J^{\pi}=(0,0)0^{+}$, interacting via various
         NN potentials. Also given are results from Ref.
         \protect\cite{GMNBPW96} with state independent TICI$_{SI}$
         and state dependent TICI $_{SD}$ correlations. }
\begin{center}
\begin{tabular}{lcccccc}
  $K_{max}$ & $N_{HH}$ & \hspace{6mm} B1       \hspace{6mm}
		       & \hspace{6mm} MT--I/III \hspace{6mm}
		       & \hspace{6mm} MTV      \hspace{6mm}
		       & \hspace{6mm} S3      \hspace{6mm}
                       & \hspace{6mm} VV   \hspace{6mm} \\ \hline
 4 &  1  & 56.71 & 52.82 & 134.29  & 31.19 & 147.42  \\
 6 &  4  & 65.39 & 59.31 & 137.72  & 38.08 & 148.70  \\
 8 & 15  & 70.03 & 60.64 & 137.80  & 42.11 & 148.49  \\
\hline
TICI$_{SI}$  &
         & 49.18 & 46.67 & 129.25  & 26.26 &         \\
\hline
TICI$_{SD}$  &
         & 61.30 & 52.67 & 130.23  & 37.30 &         \\
\end{tabular}
\end{center}
\end{table}

\begin{table}
\caption{ \label{tb:n12} \em \\
         Binding energy of twelve nucleon system $^{12}$C,
         $(L,S)J^\pi=(0,0)0^+$, interacting via various
         NN potentials. Also given are results from Refs.
         \protect\cite{GMNBPW96} (TICI),\protect\cite{Brizzi96}
         (IDEA) and \protect\cite{BBG87} (VMC).
         The TICI results as in Table \ref{tb:n8}. }

\begin{center}
\begin{tabular}{lccccc}
                  & \hspace{6mm} B1       \hspace{6mm}
	          & \hspace{6mm} MT--I/III \hspace{6mm}
	          & \hspace{6mm} MTV      \hspace{6mm}
	          & \hspace{6mm} S3      \hspace{6mm}
                  & \hspace{6mm} VV   \hspace{6mm} \\ \hline
This work               &  96.64 & 108.34 & 437.25 & 52.42 & 494.00  \\
TICI   & 103.93 & 109.04 & 429.44 & 62.99 &         \\
IDEA   &  80.1  &        &        & 44.4  &         \\
VMC    &  82.9 $\pm$ 0.2 &  &     &       &
\end{tabular}
\end{center}
\end{table}

\subsection*{Acknowledgments}
The authors would like to thank V.D. Efros, V.B. Mandelzweig and R. Guardiola
for stimulating discussions.


text:
\section{Introduction}

Finite temperature has often a very crucial influence on the
response properties of a solid. A prominent example for this
is the electrical resistivity of perfect non-magnetic metals
and ordered compounds that only take a non-zero value with a
 characteristic temperature ($T$) dependence due to thermal lattice
 vibrations. While the Holstein transport equation \cite{Hol64,Mah00} provides
 a sound basis for corresponding calculations numerical work in this
field has been done so far either on a model level or for
 simplified situations.\cite{All71,TW77,Gri76,MH83}
In practice often the Boltzmann-formalism
is adopted using the constant relaxation time ($\tau$) approximation.
 This is a very popular approach in particular when dealing with the
Seebeck effect, as in this case $\tau$ drops out.\cite{OYAA09,SY13}
The constant relaxation time approximation has also been used
extensively when dealing with the Gilbert damping
parameter $\alpha$.\cite{SF05,Kam07,GIS08} Within the description of
Kambersky \cite{Kam76,Kam07}
the conductivity- and resistivity-like
intra- and inter-band contributions to $\alpha$ show a
different dependency on $\tau$ leading
typically to a minimum for
$\alpha(\tau)$ or equivalently for $\alpha(T)$.\cite{Kam07,GIS08}
A scheme
 to deal with the temperature dependent resistivity that
is formally much more satisfying than the constant relaxation
 time  approximation is achieved by combining the
Boltzmann-formalism with a detailed calculation
of the phonon properties.
As was shown by various authors,\cite{ABK+86,SS96,XV13,XV14}
this parameter-free approach leads for non-magnetic metals in
 general to a very good agreement with experimental data.

As an alternative to this approach, thermal
lattice vibrations have also been accounted for
within various studies by quasi-static lattice
displacements leading to thermally induced structural
disorder in the system. This point of view provides the
 basis for the use of the alloy analogy, i.e.\ for
the use of
 techniques to deal with substitutional chemical disorder
also when dealing with temperature
dependent
quasi-static random  lattice displacements.
 An example for this are
 investigations on the temperature dependence of the
resistivity and the
 Gilbert parameter $\alpha$ based on the scattering
matrix approach applied to layered systems.\cite{LSYK11}
The necessary average over many configurations
of lattice displacements was taken by means of the
 super cell technique. In contrast to this the
 configurational average was determined using the
Coherent Potential Approximation (CPA) within
investigations using a Kubo-Greenwood-like linear
expression for $\alpha$.\cite{EMKK11}
 The same approach to deal with the lattice
displacements was also used recently
within calculations of angle-resolved photo emission
 spectra (ARPES) on the basis of the
one-step model of photo emission.\cite{BMM+13}

Another important contribution to the resistivity in the
case of magnetically ordered solids are thermally
induced spin fluctuations.\cite{KDT+12} Again, the alloy
analogy has been exploited extensively in the past
 when dealing with the impact of spin fluctuations
 on various response quantities. Representing a
frozen spin configuration by means of super cell
calculations has been applied for calculations of the
 Gilbert parameter for $\alpha$ \cite{LSYK11} as well as the
resistivity or conductivity, respectively.\cite{LSYK11,GBK+12,KMWB14}~
Also,
the CPA has been used for calculations of $\alpha$ \cite{MKE11}
as well as the resistivity.\cite{AD93,KDT+12} A crucial point in this
context is obviously the modeling of the temperature
dependent spin configurations. Concerning this, rather
simple models have been used,\cite{MKE11} but also quite
sophisticated schemes. Here one should mention
the transfer of data from Monte Carlo simulations
based on exchange parameters calculated in an
ab-initio way \cite{LKAG87} as well as work based on the
disordered local moment (DLM) method.\cite{AD93,GPS+85}
Although, the standard DLM does not account for transversal
spin components it nevertheless allows to represent
the paramagnetic regime with no net magnetization in
a rigorous way
Also, for the magnetically ordered
regime below the Curie-temperature it could be
demonstrated that the uncompensated
DLM (uDLM) leads for many situations still to
good agreement with experimental data on the
so-called spin disorder contribution to the resistivity.\cite{KDT+12,AD93}

In the following we present technical details and
extensions of a scheme that was already used before
 when dealing with the temperature dependence of
response quantities on the basis of Kubo's response
formalism. Various applications will be presented
for the conductivity and Gilbert
damping parameter accounting simultaneously for various types of disorder.

\section{Theoretical framework}

\subsection{Configurational average for linear response functions}

Many important quantities in spintronics can be
 formulated by making use of linear
 response formalism. Important examples for this are
the electrical conductivity,\cite{Vel69,But85}
 the spin conductivity \cite{LGK+11} or the
 Gilbert damping parameter.\cite{BTB08,EMKK11}
Restricting here for the sake of brevity to the symmetric part of the
corresponding response tensor $\chi_{\mu\nu}$
this can be expressed by a correlation function of the form:
\begin{equation}
\label{eq:sig}
\chi_{\mu\nu}
\propto \, {\rm{Tr}}\,
\big\langle \hat{A}_{\mu} \, \Im G^+ \, \hat{A}_{\nu} \, \Im G^+
\big\rangle_{\rm c} \;.
\end{equation}
It should be stressed that this not a real restriction as the
scheme described below has been used successfully when
dealing with the impact of finite temperatures
 on the anomalous Hall conductivity of Ni.\cite{KCME13}
In this case  the more complex Kubo-St\v{r}eda- or
Kubo-Bastin formulation for the
full response tensor has to be used.\cite{CB01}

The vector operator $\hat{A}_{\mu}$  in  Eq.~(\ref{eq:sig}) stands
for example in case of the electrical conductivity $\sigma_{\mu\nu}$  for the
current density operator   $\hat{j}_{\mu}$  \cite{But85}
while in case of the Gilbert damping parameter
$\alpha_{\mu\nu}$ it stands for the torque operator
 $\hat{T}_{\mu}$.\cite{SF05,EMKK11}
Within the Kubo-Greenwood-like equation
(\ref{eq:sig}) the electronic structure of the investigated
system is represented
 in terms of its retarded Green function
$ G^+ ( \vec r,\vec r\,',E) $.
 Within  multiple scattering theory
or the KKR (Korringa-Kohn-Rostoker) formalism,
$G^+(\vec{r},\vec r\,',E)$ can be written as:\cite{FS80,Wei90,ebe00}
\begin{eqnarray}
\label{eq:GFUN}
 G^+(\vec{r},\vec r\,',E)& = &
\sum_{\Lambda \Lambda'}
Z^{m}_{\Lambda}(\vec r,E)
 \tau^{mn}_{\Lambda \Lambda'}(E)
Z^{n\times}_{\Lambda'}(\vec r\,',E)\\
&&-
 \delta_{mn} \sum_{\Lambda }
   Z^{n}_{\Lambda}(\vec r,E) J^{n\times}_{\Lambda'}(\vec r'_{} ,E) \Theta(r_n'-r_n)
 \nonumber \\
&& \hspace{8.3ex}
+  J^{n}_{\Lambda}(\vec r,E) Z^{n\times}_{\Lambda'}(\vec r' ,E) \Theta(r_n-r_n')
\;.
 \nonumber
\end{eqnarray}
Here  $\vec{r},\vec{r}'$ refer  to points  within atomic  volumes around
sites $\vec R_m,  \vec R_n$, respectively, with $Z^{n}_{\Lambda}(\vec
r,E)=Z_{\Lambda}(\vec r_n,E) = Z_{\Lambda}(\vec r-\vec R_n ,E) $ being a
function centered at site $\vec R_n$.
Adopting a fully relativistic
formulation\cite{Wei90,ebe00} for Eq.~(\ref{eq:GFUN})
one gets   in a natural way  access to all
spin-orbit induced properties as for example
the anomalous and spin Hall conductivity\cite{CB01,LKE10b,LGK+11}
or Gilbert damping parameter.\cite{EMKK11}
In this case,
 the functions $Z^{n}_{\Lambda}$ and $J^{n}_{\Lambda}$ stand for the
regular and irregular, respectively,
solutions to the single-site Dirac equation
for site $n$ with the associated
single-site scattering t-matrix $ t^n_{\Lambda \Lambda'}$.
The corresponding scattering path operator
$\tau^{nn'}_{\Lambda \Lambda'}$
accounts for all scattering events connecting
 the sites $n$ and $n'$. Using a suitable spinor
 representation for the basis functions the combined quantum number
 $\Lambda=(\kappa,\mu)$  stands for the
 relativistic spin-orbit and magnetic quantum numbers $\kappa$ and $ \mu$,
respectively.\cite{Ros61,Wei90,ebe00}

As was demonstrated by various authors \cite{Vel69,But85,TKD+02b}
representing the electronic structure
 in terms of the Green function $ G^+ ( \vec r,\vec r\,',E) $
allows to account
for chemical disorder in a random
alloy by making use of a suitable alloy theory.
In this case  $\langle ... \rangle_{\rm c}$
 stands for the configurational average for a substitutional alloy
concerning the site occupation.
Corresponding expressions for the conductivity tensor
have been worked out by
 Velick\'y \cite{Vel69} and  Butler \cite{But85} using
 the single-site Coherent Potential Approximation
(CPA) that include in particular the so-called vertex corrections.

The CPA can be used to deal with chemical
but also with any other type of disorder.
In fact, making use of the different  time scales connected
with the electronic propagation and spin fluctuations the alloy analogy
is exploited when dealing with finite temperature magnetism
on the basis of the disordered local moment (DLM) model.\cite{GPS+85,SG92b}
Obviously,  the same approach can be used when dealing with
response tensors at  finite temperatures. In connection
with the conductivity
this is often called adiabatic approximation.\cite{JM90}
Following this philosophy,
the CPA has been used recently also when calculating response tensors
using Eq.~(\ref{eq:sig}) with disorder in the system caused by
thermal lattice vibrations \cite{EMKK11,KCME13}
as well as spin fluctuations.\cite{KDT+12,MKWE13}

\subsection{Treatment of thermal lattice displacement}
\label{sec-thermal-lattice-displacements}
A way to account for the  impact of the thermal displacement of atoms
from their equilibrium positions,
i.e.\ for thermal lattice vibrations,
on the electronic structure is to set up
a representative  displacement configuration for the
atoms within an enlarged unit cell (super-cell  technique).
In this case one has to use either a
very large super-cell or to take the average
over a set of super-cells.
Alternatively, one may make use of the alloy analogy
for the averaging problem.
This allows in particular to restrict
to the standard unit cell.
Neglecting the correlation
 between the thermal displacements of neighboring
atoms from their equilibrium positions
the properties of the thermal averaged system
can be deduced by making use of the single-site
CPA. This basic idea is illustrated by
Fig.\ \ref{CPA-displacement}.
\begin{figure}[hbt]
   \includegraphics[angle=0,width=0.9\linewidth,clip]{Fig1_vib_A4.eps}
\caption{\label{CPA-displacement}
Configurational averaging for thermal lattice displacements:
the continuous distribution
$P(\Delta \vec{R}_n(T))$  for the atomic displacement
vectors is replaced by a discrete set of vectors
$\Delta \vec{R}_v(T)$ occurring
with the probability $x_v$. The configurational average for this
 discrete set of displacements is made using the CPA leading
to a periodic effective medium.}
\end{figure}
To make use of this scheme a discrete
set of $N_{v}$ displacement vectors
$\Delta \vec{R}^q_v(T)$
with probability $x^q_v$
($v=1,..,N_{v}$)
is constructed   for each basis atom $q$
within the standard unit cell that is conform with the local symmetry and
 the temperature dependent root mean  square
displacement $(\langle u^2\rangle_T)^{1/2}$ according to:
\begin{equation}
\label{eq:displacement}
\frac{1}{N_{v}}
\sum_{v=1}^{N_{v}} | \Delta \vec{R}^q_v(T) |^2 = \langle u_q^2\rangle_T \;.
\end{equation}
In the general case, the mean square displacement along the direction
  $\mu$ ($\mu = {x,y,z}$) of the atom $i$ can be either taken from
experimental data or represented by the expression based on the
phonon calculations \cite{Boe83}

\begin{equation}
\label{eq:MSD}
\langle
 u^2_{i,\mu}
\rangle_T = \frac{3\hbar}{2M_i}
\int_0^{\infty} d\omega
g_{i,\mu}(\omega)\frac{1}{\omega}
{\rm{coth}}\frac{\hbar\omega}{2k_{\rm B}T}  \;,
\end{equation}
where
$h=2\pi\hbar$ the Planck constant,  $k_{\rm B}$ the Boltzmann constant,
$g_{i,\mu}(\omega)$ is a partial phonon density of states. \cite{Boe83}
On the other hand, a rather good estimate for the root mean square
displacement can be obtained using  Debye's theory. In this
case, for systems with one atom per unit cell, Eq.\ (\ref{eq:MSD}) can
be reduced to the expression:
\begin{equation}
\label{eq:Debye}
\langle u^2\rangle_T =
\frac{1}{4}\frac{3h^2}{\pi^2Mk_{\rm B}
\Theta_D}\left[\frac{\Phi(\Theta_D/T)}{\Theta_D/T}
+ \frac{1}{4}\right]
\end{equation}
with $\Phi(\Theta_D/T)$ the Debye function and
$\Theta_D$ the Debye temperature \cite{GMMP83}.
Ignoring the zero temperature term $1/4$ and
assuming a frozen potential for the atoms, the situation can be dealt
with in full analogy to the treatment of disordered alloys on the basis
of the CPA. The probability $x_v$ for a specific displacement
$v$ may normally be chosen as $1/N_{v}$.
The Debye temperature $\Theta_D$ used in Eq.~(\ref{eq:Debye}) can be
either taken from  experimental data or calculated
representing it in terms of the elastic constants \cite{FBS01}. In general
the latter approach
should give more reliable results in the case of multicomponent systems.

To simplify notation we restrict in the following to
systems with one atom per unit cell.
The index $q$ numbering sites in the unit cell can therefore
be dropped, while the  index $n$ numbers the lattice sites.

Assuming a rigid displacement of the atomic potential
in the spirit of the rigid muffin-tin approximation \cite{PZDS97,Lod76}
the corresponding single-site t-matrix $ \underline{t}^{\rm loc}$
with respect to the local frame of reference
connected with the displaced atomic position
is unchanged.
With respect to the global frame of  reference
connected with the equilibrium atomic positions $ \vec{R}_n$, however,
the corresponding t-matrix $ \underline{t}$
is given by the transformation:
\begin{equation}
\label{eq:U-trans}
\underline{t} = \underline{U}(\Delta \vec{R})\,\underline{t}^{\rm loc}\,
                     \underline{U}(\Delta \vec{R})^{-1} \;.
\end{equation}
The so-called U-transformation matrix $\underline{U}(\vec{s})$
is given  in its
non-relativistic form
by:\cite{Lod76,PZDS97}
\begin{equation}
\label{eq:U-trans-matrix}
U_{LL'}(\vec{s}) =
4\pi \sum_{L''}i^{l+l''-l'}\, C_{LL'L''}\, j_{l''}(|\vec{s}|k)\, Y_{L''}(\hat{s})
\;.
\end{equation}
Here $L=(l,m)$ represents the non-relativistic
angular momentum quantum numbers,
$j_{l}(x)$ is a spherical Bessel function,
$Y_{L}(\hat{r})$ a real spherical harmonics,
$C_{LL'L''}$ a corresponding  Gaunt number
and $k=\sqrt{E}$ is the electronic wave vector.
The relativistic version of the U-matrix
is obtained by a standard Clebsch-Gordan transformation.\cite{Ros61}

The various displacement vectors $\Delta \vec{R}_v(T)$
can be used to determine the properties of a
 pseudo-component of a pseudo alloy.
Each of the
 $N_{v}$ pseudo-components with $|\Delta \vec{R}_v(T)| = \langle
 u^2\rangle_T ^{1/2}$ is characterized
by a corresponding  U-matrix
 $\underline{U}_v$ and
t-matrix  $\underline{t}_v$. As for a substitutional alloy
the configurational
average can be determined by solving the
multi-component CPA equations
within the global frame of  reference:
\begin{eqnarray}
\label{eq:CPA1}
\underline{\tau}_{{\rm CPA}}^{nn} &= &
\sum_{v=1}^{N_{v}}
x_v \underline{\tau}_{v}^{nn}
\\
\label{eq:CPA2}
\underline{\tau}_{v}^{nn}& = &
\big[
    (\underline{t}_{v})^{-1}
-   (\underline{t}_{{\rm CPA}})^{-1}
+   (\underline{\tau}_{{\rm CPA}}^{nn})^{-1}
\big]^{-1}
\\
\label{eq:CPA3}
\underline{\tau}_{{\rm CPA}}^{nn}
 & = & \frac{1}{\Omega_{{\rm BZ}}}
 \int_{\Omega_{\rm BZ} } d^{3}k
\left[ (\underline{t}_{{\rm CPA}})^{-1}
      - \underline{G}(\vec{k},E)  \right]^{-1}  \; ,
\end{eqnarray}
where the underline indicates
matrices with respect to the
combined index $ \Lambda$.
As it was pointed out in the previous work \cite{MKWE13},
the cutoff for the angular momentum expansion in these
calculations should be taken $l \geq l_{max} + 1$ with the $l_{max}$
value used in the calculations for the non-distorted lattice.

The first of these CPA equations represents the
 requirement for the mean-field CPA medium
that  embedding of a component $v$
should lead in the average to no additional scattering.
Eq.~(\ref{eq:CPA2}) gives the scattering path operator for the
embedding of the component $v$ into the
CPA medium while Eq.~(\ref{eq:CPA3})
gives the  CPA scattering path operator
in terms of a Brillouin
zone integral with $ \underline{G}(\vec{k},E)$ the so-called
KKR structure constants.

Having solved the CPA equations
the linear response quantity of interest
may be calculated using Eq.~(\ref{eq:sig}) as for an
ordinary substitutional alloy.\cite{Vel69,But85}
This implies that one also have to
deal with the so-called vertex corrections\cite{Vel69,But85}
that take into account that one has to
deal with a configuration average of the type
$\langle \hat{A}_{\mu} \, \Im G^+ \,
         \hat{A}_{\nu} \, \Im G^+
 \rangle_{\rm c} $
 that in general will differ from the
simpler product
$\langle \hat{A}_{\mu} \, \Im G^+ \,
 \rangle_{\rm c}
 \langle
          \hat{A}_{\nu} \, \Im G^+
 \rangle_{\rm c} $.

\subsection{Treatment of thermal spin fluctuations}

As for the disorder connected with
thermal displacements the impact of disorder
due to thermal spin fluctuations
may be accounted for by use of the super-cell
technique.~
Alternatively one may again use the alloy analogy
and determine the necessary configurational average
by means of the CPA
as indicated
in Fig.\ \ref{CPA-fluctuation}.
\begin{figure}[hbt]
 \begin{center}
   \includegraphics[angle=0,width=0.9\linewidth,clip]{Fig2_fluct_A4.eps}
 \end{center}
  \caption{
\label{CPA-fluctuation}
Configurational averaging for thermal spin fluctuations:
the continuous distribution $P(\hat{e}_n)$ for the orientation of
 the magnetic moments is replaced by a discrete set of
orientation vectors $\hat{e}_f$ occurring with a probability $x_f$.
 The configurational average for this discrete set of
orientations is made using the CPA leading to a periodic effective medium.
}
\end{figure}
As for the thermal displacements
in a first step a set of representative
orientation vectors $\hat{e}_f$
(with $f=1,...,N_f$) for the local magnetic
moment is introduced (see below).
Using the rigid spin approximation
the spin-dependent
part  $B_{\rm xc}$ of the exchange-correlation potential
does not change for the local frame of  reference fixed to
the magnetic moment
when the moment is oriented along an
orientation vector  $\hat{e}_f$ .
This implies that the single-site t-matrix
 $ \underline{t}_f^{\rm loc}$ in the local frame is the same
for all orientation vectors. With respect
to the common global frame
that is used to deal with the multiple scattering
(see Eq.~(\ref{eq:CPA3}))
the t-matrix for a given orientation vector
is determined by:
\begin{equation}
\label{eq:R-trans}
\underline{t} = \underline{R}(\hat{e})\,\underline{t}^{\rm loc}\,
                \underline{R}(\hat{e})^{-1} \;.
\end{equation}
Here the transformation from the
local to the global frame of reference
is expressed by the rotation matrices $ \underline{R}(\hat{e})$
that are determined by the vectors $\hat{e}$
or corresponding Euler angles.\cite{Ros61}

Again the configurational average
for the pseudo-alloy can be obtained
by setting  up and solving
CPA equations in analogy to
Eqs.~(\ref{eq:CPA1}) to (\ref{eq:CPA3}).

\subsection{Models of spin disorder}
\label{seq:Model-spin-disorder}

The central problem with the scheme described above is obviously
 to construct a realistic and
 representative set of orientation vectors $\hat{e}_f$
and probabilities  ${x}_f$
 for each temperature $T$. A rather appealing approach is to calculate
the exchange-coupling parameters $J_{ij}$ of a system
 in an ab-initio way \cite{LKAG87,USPW03,EM09a}
and to use them in subsequent Monte Carlo simulations.
Fig.\ \ref{FIG-MAG-T} (top) shows results for
the temperature dependent  average reduced
magnetic moment
of  corresponding simulations for bcc-Fe
 obtained for a periodic cell with 4096 atom sites.
\begin{figure}[hbt]
 \begin{center}
   \includegraphics[angle=0,width=0.7\linewidth,clip]{Fig_3a_Fe_M_T_atoms4096_5a.eps}
   \includegraphics[angle=0,width=0.7\linewidth,clip]{Fig_3b_Fe_M_T_fit_to_MC_mod.eps}
   \includegraphics[angle=0,width=0.7\linewidth,clip]{Fig_3c_Fe_M_T_smooth_5b.eps}
  \caption{\label{FIG-MAG-T}
Averaged reduced magnetic
moment
$M(T) = \langle m_z \rangle_T / |\langle \vec{m} \rangle_{T=0}|$
along the z-axis as a function of the
temperature $T$.
Top: results of Monte Carlo simulations using
scheme MC* (full squares) compared with results of
subsequent KKR-calculations (open squares).
Middle: results of Monte Carlo simulations using
scheme MC (full squares) compared with results
using a mean-field fit with a constant Weiss field
$w_{\rm MC}(T_C)$ (open diamonds) and a temperature dependent
Weiss field $w_{\rm MC}(T)$ (open squares).
In addition experimental data (full circles) together
with a corresponding  mean-field fit obtained for a
temperature dependent
Weiss field $w_{\rm exp}(T)$.
Bottom: results of Monte Carlo simulations using
scheme MC (full squares) compared with results
subsequent KKR-calculations
using the MC (triangles up)
and a corresponding DLM (triangle down) spin configuration,
respectively.
}
 \end{center}
\end{figure}
The full line gives the
value for the reduced magnetic moment
$M_{\rm MC^*}(T) = \langle m_z \rangle_T /m_0$
projected on the z-axis for the last Monte Carlo step
($\hat z$ is the orientation of the total moment, i.e.\
$\langle \vec{m}\rangle_T \| \hat z $; the saturated magnetic moment at
$T = 0$~K is $m_0 = |\langle \vec{m} \rangle_{T=0}|$). This scheme is
 called ${\rm MC^*}$ in the following.
In spite of the rather large number of sites (4096) the curve
 is rather noisy in particular when
approaching the Curie temperature.
Nevertheless, the spin configuration of
the last MC step was used as an input for subsequent
SPR-KKR-CPA calculations using
 the orientation vectors $\hat{e}_f$ with the probability $x_f=1/N_f$
 with $N_f = 4096$.
As Fig.\ \ref{FIG-MAG-T} (top) shows, the
temperature dependent reduced magnetic moment $M_{\rm KKR(MC^*)}(T)$ deduced
from the electronic structure calculations follows one-to-one the
 Monte Carlo data $M_{\rm MC^*}(T)$. This is a very encouraging result for
further applications (see below) as it demonstrates that the
CPA although being a mean-field method and used here in its single-site
formulation is nevertheless capable to reproduce results of MC simulations
 that go well beyond  the mean-field level.

However, using the set of vectors $\hat{e}_f$ of scheme
MC* also for calculations of the
Gilbert damping parameters $\alpha$ as a
function of temperature led to extremely
 noisy and unreliable curves for $\alpha(T)$. For that reason an average has
been taken over many MC steps
(scheme MC) leading to a much smoother
 curve for $M_{\rm MC}(T)$ as can be seen from
Fig.\ \ref{FIG-MAG-T} (middle) with
 a Curie temperature $T_{\rm C}^{\rm MC}=1082$~K.
 As this enlarged set of vectors
$\hat{e}_f$ got too large to be used directly  in subsequent
SPR-KKR-CPA calculations, a scheme was worked out
 to get a set of vectors $\hat{e}_f$ and
probabilities $x_f$ that is not too large
but nevertheless leads to smooth curves
for $M(T)$.

The first attempt was to use
 the Curie temperature $T_{\rm C}^{\rm MC}$ to
deduce a corresponding temperature independent Weiss-field $w(T_{\rm C})$
on the basis of the standard
 mean-field relation:
\begin{equation}
\label{eq:wTC}
w(T_{\rm C}) = \frac{3k_{\rm B}T_C}{m_0^2} \;.
\end{equation}
 This leads to a reduced magnetic moment
curve $M_{\rm MF}(T)$ that shows by construction the same Curie temperature as
 the MC simulations.
For temperatures between $T=0$~K  and $T_{\rm C}$, however, the
 mean-field reduced magnetic moment $M_{\rm MF}(T)$ is well below the MC curve
(see Fig.\ \ref{FIG-MAG-T} (middle) ).

As an
alternative to this simple approach we introduced a temperature dependent
 Weiss field $w(T)$.
This allows to describe the temperature dependent magnetic properties
using the results obtained beyond the mean-field approximation.
At the same time the calculation of the statistical average can be performed
treating the model Hamiltonian in terms of the mean field theory.
For this reason the reduced magnetic moment  $M(T)$, being a solution
of equation (see e.g.\ \cite{Tik64})
\begin{equation}
\label{eq:xf}
 M(T)  = L\left(\frac{w m_0^2 M(T)}{k_{\rm B}T}\right) \;,
\end{equation}
was fitted to that obtained  from MC simulations $M_{\rm MC}(T)$
with the Weiss field $w(T)$ as a fitting parameter, such that
\begin{equation}
\label{eq:xf-a}
\lim_{w \to w(T)} M(T)  = M_{MC}(T) \;,
\end{equation}
with $L(x)$ the Langevin function.

The corresponding temperature dependent probability
 $x(\hat{e})$ for an atomic magnetic moment to be oriented along $\hat{e}$ is
proportional to $\exp(-w(T) \hat z \cdot \hat e/k_{\rm B}T)$ (see,
e.g. \cite{Tik64}).
 To calculate this value we
used $N_{\theta}$ and $N_{\phi}$ points for
 a regular grid for the spherical angles $\theta$ and $\phi$ corresponding
to the vector $\hat{e}_f$:
\begin{equation}
\label{eq:xf-b}
       x_f = \frac{\exp(-w(T) \hat z \cdot \hat{e}_{f}/k_{\rm B}T)}
  {\sum_{f'} \exp(-w(T) \hat z \cdot \hat{e}_{f'}/k_{\rm B}T)}   \;.
\end{equation}
 Fig.\ \ref{FIG-MAG-PROB} shows for three different
temperatures the $\theta$-dependent behavior of $x(\hat{e})$.
\begin{figure}[hbt]
 \begin{center}
   \includegraphics[angle=0,width=0.7\linewidth,clip]{Fig_4a_Fe_M_THETA_200K.eps}
   \includegraphics[angle=0,width=0.7\linewidth,clip]{Fig_4b_Fe_M_THETA_400K.eps}
   \includegraphics[angle=0,width=0.7\linewidth,clip]{Fig_4c_Fe_M_THETA_800K.eps}
 \caption{\label{FIG-MAG-PROB}
Angular distribution $P(\theta) $
of the atomic magnetic moment $\vec{m}$
obtained from Monte Carlo simulations (MC)
for the temperature $T= 200 $, $400 $, and $800$~K
compared with field mean-field (MF) data, $x_f$, (full line)
obtained by fitting using a
temperature dependent Weiss field $w(T)$ (Eq.\ \ref{eq:xf}).
}
 \end{center}
\end{figure}
As one notes, the MF-fit to the MC-results perfectly reproduces
 these data for all temperatures.
 This applies of course not only for the angular resolved
 distribution of the magnetic moments shown in
Fig.\ \ref{FIG-MAG-PROB} but also for the average
reduced magnetic moment recalculated using Eq.(\ref{eq:xf}), shown in
Fig.\ \ref{FIG-MAG-T}. Obviously, the MF-curve $M_{\rm MF(MC)}(T)$
obtained using the temperature dependent Weiss field parameter $w(T)$ perfectly reproduces the
 original $M_{\rm MC}(T)$ curve. The great advantage
of this fitting procedure is that it
allows to replace the MC data set with a large
 number $N_f^{\rm MC}$ of orientation vectors
$\hat{e}_f$ (pointing in principle into any direction)
 with equal probability $x_f=1/N_f^{\rm MC}$
by a much smaller data set with
$N_f= N_\theta \,  N_\phi$
with $x_f$
 given by Eq.~(\ref{eq:xf-b}).

Accordingly, the reduced
 data set can straight forwardly
 be used for subsequent electronic structure calculations.
Fig.\ \ref{FIG-MAG-T} (bottom) shows that the calculated
temperature dependent
 reduced magnetic moment $M_{\rm KKR-MF(MC)}(T)$ agrees perfectly
 with the reduced magnetic moment
$M_{\rm MC}(T)$ given by the underlying MC simulations.

The DLM method has the appealing feature that it combines
 ab-initio calculations
 and thermodynamics in a coherent way.
 Using a non-relativistic formulation, it
was shown that the corresponding averaging
over all orientations of the individual
atomic reduced magnetic moments can be mapped onto a
 binary pseudo-alloy with one
pseudo-component having up- and
 downward orientation of the spin moment  with concentrations
$x_{\uparrow}$ and $x_{\downarrow}$, respectively.\cite{AD93,Aka98}
 For a fully relativistic formulation,
with spin-orbit coupling included, this
simplification cannot be justified anymore
 and a proper average has to be taken over
 all orientations.\cite{SOR+04}
As we do not perform DLM calculations but use
here only the DLM picture to
represent MC data, this complication is
ignored in the following. Having the set of
orientation vectors $\hat{e}_f$ determined by
 MC simulations the corresponding
concentrations $x_{\uparrow}$ and $x_{\downarrow}$
 can straight forwardly be fixed for each
temperature by the requirement:
\begin{equation}
\label{eq:DLM-requ}
\frac{1}{N_f} \sum_{f=1}^{N_f}
\hat{e}_f =
x_{\uparrow} \hat z +  x_{\downarrow}(- \hat z )
 \;,
\end{equation}
with $x_{\uparrow} + x_{\downarrow} = 1$.
Using this simple scheme electronic structure calculations have been performed
for a binary alloy having collinear magnetization. The resulting
reduced magnetic moment $M_{\rm KKR-DLM(MC)}(T)$ is shown in
Fig.\ \ref{FIG-MAG-T} (bottom). As one notes, again the original MC
results are perfectly reproduced. This implies that when calculating
the projected reduced magnetic moment $M_{\rm z}$ that is determined by
the averaged Green function $\langle G \rangle$  the
transversal magnetization has hardly any impact.

Fig.\ \ref{FIG-MAG-T} (middle) gives also
experimental data for the $M(T)$.\cite{CG71}
While the experimental Curie-temperature
 $T_{\rm C}^{\rm exp}=1044$~K\cite{CG71}
is rather well reproduced by the
MC simulations $T_{\rm C}^{\rm MC}=1082$~K
one notes that the MC-curve $M_{\rm MC}(T)$
is well below the experimental curve. In particular,
$M_{\rm MC}(T)$ drops too
fast with increasing $T$ in the low temperature regime and does not show the
$T^{3/2}$-behavior.
The reason for this is that the MC simulations do not
properly account for the low-energy long-ranged spin wave excitations
responsible for the low-temperature magnetization variation.
 Performing ab-initio calculations for the spin wave energies and using
 these data for the calculation of $M(T)$ much better agreement with
experiment can indeed be obtained in the low-temperature regime than with
 MC simulations.\cite{HEPO98}

As the fitting scheme sketched above needs only
the temperature reduced magnetic moment $M(T)$ as input it can be applied not only
to MC data but also to experimental data.
Fig.\ \ref{FIG-MAG-T} shows that the mean field fit $M_{\rm MF(exp)}(T)$ again
perfectly fits the experimental
 reduced magnetic moment curve $M_{\rm exp}(T)$. Based on
this good agreement this
corresponding data set $\{\hat{e}_f,x_f\}$
has also been used
for the calculation of response tensors (see below).

An additional much
simpler scheme to simulate the experimental $M_{\rm exp}(T)$ curve is to assume the
individual atomic moments to be distributed on a cone, i.e.\  with $N_\theta=1$
 and $N_{\phi}>>1$.\cite{MKE11}
 In this case the opening angle $\theta(T)$ of the cone is chosen
such to reproduce $M(T)$. In contrast to the standard DLM picture,
this simple scheme allows already to account for transversal components
 of the magnetization. Corresponding results for response tensor calculations will be shown below.

Finally, it should be stressed here that the various spin configuration models
discussed above assume a rigid spin moment, i.e.\ its magnitude does not change with
temperature nor with orientation. In contrast to this
Ruban et al.\ \cite{RKMJ07}
use a longitudinal spin fluctuation Hamiltonian with the corresponding parameters
derived from ab-initio calculations.
As a consequence, subsequent Monte Carlo simulations based on this Hamiltonian
account in particular for longitudinal fluctuations of the spin moments.
A similar approach has been used by Drchal et al. \cite{DKT13,DKT13a}
leading to good agreement with the results of Ruban et al.
However, the scheme used in these calculations does not supply in a
straightforward manner the necessary input for temperature dependent
transport calculations.
This is different from the work of Staunton et al. \cite{SBS+14} who
performed self-consistent relativistic DLM calculations without the
restriction to a collinear spin configuration. This approach in
particular accounts in a self-consistent way for longitudinal spin
fluctuations.

\subsection{Combined chemical and thermally induced disorder}

The various types of disorder discussed above may be combined
with each other as well as with chemical i.e.\  substitution
disorder. In the most general case a pseudo-component $(vft)$
is characterized by its chemical  atomic type $t$, the spin
fluctuation $f$ and lattice displacement $v$. Using the rigid
muffin-tin and rigid spin approximations, the single-site
t-matrix $t^{\rm loc}_t $ in the local frame is independent from the
orientation vector $\hat{e}_f$ and displacement vector
$\Delta \vec{R}_{v}$,
and coincides with $ \underline{t}_t$ for the atomic type $t$. With
respect to the common global frame one has accordingly
 the t-matrix:
\begin{equation}
\label{eq:tvft}
 \underline{t}_{vft} =     \underline{U}(\Delta \vec{R}_{v})\,
                     \underline{R}(\hat{e}_{f})\,
                     \underline{t}_t \,
                     \underline{R}(\hat{e}_{f})^{-1}
                     \underline{U}(\Delta \vec{R}_{v})^{-1}
\;.
\end{equation}
With this the corresponding CPA equations are identical
 to
Eqs.~(\ref{eq:CPA1}) to (\ref{eq:CPA3})
 with the index $v$ replaced by the
combined index $(vft)$. The corresponding
 pseudo-concentration $x_{vft}$ combines the concentration
 $x_t$ of the atomic type $t$ with the probability for
the orientation vector $\hat{e}_f$ and displacement vector $\Delta \vec{R}_{v}$.

\section{Computational details}

The electronic structure of the investigated ferromagnets
bcc-Fe   and fcc-Ni,
has been calculated self-consistently using the
spin-polarized  relativistic
KKR (SPR-KKR)
band structure
 method.\cite{EKM11,SPR-KKR6.3}
For the exchange correlation potential the parametrization
as given by Vosko et al.\ \cite{VWN80} has been used.
The angular-momentum cutoff of $l_{max} = 3$  was used in the KKR multiple
scattering expansion.
The lattice parameters
 have been set to the experimental values.

In a second step the
 exchange-coupling parameters $J_{ij}$ have been calculated
 using the so-called Lichtenstein formula.\cite{LKAG87} Although the
SCF-calculations have been done on a fully-relativistic level the
 anisotropy of the exchange coupling due to the spin-orbit coupling has
 been neglected here. Also, the small influence of the magneto-crystalline
 anisotropy for the subsequent Monte Carlo (MC) simulations has been
ignored, i.e.\  these have been based on a classical Heisenberg Hamiltonian.
The MC simulations were done in a standard way using the Metropolis algorithm
 and periodic boundary conditions.
 The theoretical Curie temperature $T_{\rm C}^{\rm MC}$ has been
deduced from the maximum of the magnetic susceptibility.

The temperature dependent spin configuration obtained during a MC
simulation has been used to construct
a set of orientations $\hat{e}_f$
and probabilities  ${x}_f$ according to the  schemes
MC* and MC described in section \ref{seq:Model-spin-disorder}
to
be used within subsequent SPR-KKR-CPA calculations
 (see above). For the corresponding
calculation of the reduced magnetic moment the potential obtained from the
SCF-calculation for the perfect ferromagnetic state ($T=0$K) has been used.
The calculation for the electrical conductivity as well as the Gilbert
 damping parameter has been performed as
described elsewhere.\cite{KLSE11,MKWE13}

\section{Results and discussion}

\subsection{Temperature dependent conductivity}

  Eq.~(\ref{eq:sig})  has been used together with the various
schemes described above to calculate the temperature
dependent longitudinal resistivity  $ \rho (T) $
 of the pure
ferromagnets Fe, Co and Ni. In this case obviously
disorder due to thermal displacements of  the atoms
as well as spin fluctuations contribute to the resistivity.

To give an impression on the impact of the
thermal displacements alone
  Fig.\ \ref{FIG-RHO-Cu}
gives the temperature dependent
 resistivity  $ \rho (T) $
of  pure Cu ($ \Theta_{\rm Debye}=315$ K)
that is found in very good agreement with
corresponding experimental data.\cite{Bas82}
\begin{figure}[hbt]
 \begin{center} \includegraphics[angle=0,width=0.8\linewidth,clip]{Fig_5_Cu_sigma_T.eps}
  \caption{\label{FIG-RHO-Cu}
Temperature dependent longitudinal
 resistivity of fcc-Cu $ \rho(T)$
obtained by accounted for thermal vibrations
as described in section \ref{sec-thermal-lattice-displacements} compared
with corresponding experimental data.\cite{Bas82}
 In addition results are shown based on the LOVA (lowest order variational approximation) to the
Boltzmann formalism.\cite{SS96}
}
 \end{center}
\end{figure}
This implies that the alloy analogy model
that ignores any inelastic scattering events
should in general lead to rather reliable results for the
resistivity induced by thermal displacements.
Accordingly, comparison with experiment should allow
for magnetically ordered systems
to find out the most appropriate model
for spin fluctuations.

 Fig.\ \ref{FIG-RHO-DISP-FLUCT} (top)
shows theoretical results for $ \rho (T) $ of bcc-Fe
due to thermal displacements $ \rho_{v}(T) $,
spin fluctuations described by the scheme MC
$ \rho_{ \rm MC}(T)$ as
 well as the combination of the two
influences ($ \rho_{ \rm v,MC}(T)$).
\begin{figure}[hbt]
 \begin{center}
   \includegraphics[angle=0,width=0.8\linewidth,clip]{Fig_6a_Fe_rho_xx_vib_fluct_MC.eps}
   \includegraphics[angle=0,width=0.8\linewidth,clip]{Fig_6b_Fe_fluct_MC60_DLM_vs_all_dir.eps}
  \caption{\label{FIG-RHO-DISP-FLUCT}
Temperature dependent longitudinal
 resistivity of bcc-Fe $ \rho(T)$
obtained by accounted for thermal vibrations
and spin fluctuations
as described in section \ref{sec-thermal-lattice-displacements}.
Top:  accounting for
vibrations                        (vib, diamonds),
spin fluctuations using scheme MC (fluct, squares)
and both                          (vib+fluct, circles).
Bottom: accounting for spin fluctuations $\hat{e}_f = \hat{e}(\theta_f,\phi_f)$
using the schemes:  MC       (squares) with $0 \leq \theta_f \leq \pi; 0 \leq \phi_f \leq 2\pi$,
                    DLM(MC)  (triangles up) with $\theta_{f1} = 0, \theta_{f2} = \pi$,
and                 cone(MC) (triangles down) $\theta_f = \langle \theta_f
\rangle_T; 0 \leq \phi_f \leq 2\pi$.
The full and open
symbols represent the results obtained
with the vertex corrections included (VC)
 and excluded (NV), respectively.
}
 \end{center}
\end{figure}
First of all one notes that $ \rho_{ \rm v}(T) $
is not influenced within the adopted model
by the Curie temperature $T_{ \rm C}$
 but is determined
only by the Debye temperature.
$ \rho_{ \rm MC}(T) $, on the other hand,
reaches saturation for   $T_{ \rm C}$
as the spin disorder does not increase anymore
with increasing temperature in the
paramagnetic regime.
 Fig.\ \ref{FIG-RHO-DISP-FLUCT} also
shows that  $ \rho_{ \rm v}(T) $
and
$ \rho_{ \rm MC}(T) $
are comparable for low temperatures
but $ \rho_{ \rm MC}(T) $ exceeds
 $ \rho_{ \rm v}(T) $ more and more
for higher temperatures.
Most interestingly, however,
the resistivity for the combined influence
of thermal displacements and spin fluctuations
$ \rho_{ \rm v,MC}(T) $ does
not coincide  with the sum of $ \rho_{ \rm v}(T) $
and
$ \rho_{ \rm MC}(T) $
but exceeds the sum for low temperatures
and lies below the sum when approaching
  $T_{ \rm C}$.

Fig.\ \ref{FIG-RHO-DISP-FLUCT} (bottom)
shows the results   of three  different
calculations including the effect of
spin fluctuations
as a function of the temperature.
The curve $ \rho_{ \rm MC}(T) $
is identical with that given in
Fig.\ \ref{FIG-RHO-DISP-FLUCT} (top)
based on Monte Carlo simulations.
The curves $ \rho_{ \rm DLM(MC)}(T)$
and
$ \rho_{ \rm cone(MC)}(T)$
are based on a DLM- and cone-like representation
of the MC-results, respectively.
For all three cases results are given including
as well as ignoring the vertex corrections.
As one notes the vertex corrections
play a negligible role for all three
spin disorder models.
This is fully in line with the experience for
the longitudinal resistivity of disordered
transition metal alloys: as long as the
the states at the Fermi level have
dominantly d-character the vertex corrections
can be neglected in general. On the other
hand, if the sp-character dominates inclusion of
vertex corrections may
alter the result in the order of 10 \%.\cite{BEWV94,TKDW04}

Comparing
the DLM-result $ \rho_{ \rm DLM(MC)}(T)$ with
$ \rho_{MC}(T) $
one notes
in contrast to the results for $ M(T) $
shown above (see Fig.\ \ref{FIG-MAG-T} (bottom))
quite an appreciable deviation.
This implies that the restricted
collinear representation
of the spin configuration implied by the DLM-model
introduces errors for the
configurational average that seem in general
to be unacceptable,
For the Curie temperature and beyond in the
paramagnetic regime
$ \rho_{ \rm DLM(MC)}(T)$ and
$ \rho_{ \rm MC}(T) $
coincide, as it was shown formally before.\cite{KDT+12}

Comparing finally
$ \rho_{ \rm cone(MC)}(T)$
based on the conical representation of
the MC spin configuration with
$ \rho_{ \rm MC}(T) $
one notes that also this simplification
leads to quite strong deviations from the
more reliable result.
Nevertheless, one notes that
$ \rho_{ \rm DLM(MC)}(T)$
agrees  with
$ \rho_{ \rm MC}(T) $
for the Curie temperature and also
accounts to some
extent for the impact of the transversal
components of the magnetization.

The theoretical results for bcc-Fe ($ \Theta_{\rm Debye}=420$ K) based
on the combined inclusion of the effects of
thermal displacements
and
spin fluctuations using the MC scheme
($ \rho_{ \rm v,MC}(T) $)
are compared in
Fig.\ \ref{FIG-RHO-EXPT} (top)
with experimental data ($ \rho_{ \rm exp}(T) $).
\begin{figure}[hbt]
 \begin{center}
  \includegraphics[angle=0,width=0.8\linewidth,clip]{Fig_7a_Fe_rho_xx_vib_fluct_MC_exp.eps}
  \includegraphics[angle=0,width=0.8\linewidth,clip]{Fig_7b_rho_xx_Ni_vib_fluct.eps}
 \end{center}
  \caption{\label{FIG-RHO-EXPT}
Top: Temperature dependent longitudinal
 resistivity of bcc-Fe $ \rho(T)$
obtained by accounted for thermal vibrations
and spin fluctuations
 using the scheme MC
(vib+fluct(MC), squares)
and a mean-field fit to the experimental
temperature
magnetic
moment $M_{\rm exp} $
(vib+fluct(exp), diamonds)
compared with experimental data (circles).\cite{Bas82}
Bottom: corresponding results for fcc-Ni.
In addition results are shown accounting for thermal displacements
(vib) only for the ferromagnetic (FM) as well paramagnetic (PM) regime.
Experimental data have been taken from Ref.\ \onlinecite{HAW+83}.
}
\end{figure}
For the Curie temperature obviously a very good
agreement with experiment is found
while for lower temperatures $ \rho_{ \rm v,MC}(T) $
exceeds $ \rho_{ \rm exp}(T) $.
This behavior correlates well with
that of
the temperature dependent reduced magnetic moment $M(T)$
shown in  Fig.\ \ref{FIG-MAG-T} (middle).
The too rapid decrease of $M_{ \rm MC}(T) $
compared with experiment implies
an essentially overestimated
spin disorder at any temperature leading in turn to a
too large resistivity
$ \rho_{ \rm v,MC}(T) $.
On the other hand, using the
temperature dependence of the
experimental reduced magnetic moment $M_{ \rm exp}(T) $
to set up the temperature dependent spin
configuration as described above
a very satisfying agreement is found with the
experimental resistivity data $ \rho_{ \rm exp}(T) $.
Note also that above $T_C$ the calculated resistivity
riches the saturation in contrast to the experimental data
where the continuing increase of $ \rho_{ \rm exp}(T) $ can be
attributed to the longitudinal spin fluctuations leading
to a temperature dependent distribution of local magnetic moments
on Fe atoms.\cite{RKMJ07}
However, this contribution was not taken into account
because of restriction in present calculations
using fixed value for the local reduced magnetic moments.

Fig.\ \ref{FIG-RHO-EXPT} (bottom) shows corresponding results
for the temperature dependent resistivity
of fcc-Ni ($ \Theta_{\rm Debye}=375$ K). For the ferromagnetic (FM) regime that the theoretical results
are comparable in magnitude when only thermal displacements  ($ \rho_{ \rm v}(T) $)
or spin fluctuations ($ \rho_{ \rm MF}(T) $)
are accounted for.
In the later case the mean field $w(T)$  has been fitted to the experimental
 $M(T)$-curve.    Taking both into account leads to
a resistivity  ($ \rho_{ \rm v,MF}(T) $)
 that are  well above the sum of the individual terms $ \rho_{ \rm v}(T) $ and
 $ \rho_{ \rm MF}(T) $.
Comparing $ \rho_{ \rm v,MF}(T) $ with experimental data
 $ \rho_{ \rm exp}(T) $ our finding shows that the theoretical results
overshoots the experimental one the closer one comes to the critical temperature.
This is a clear indication that the assumption of a rigid spin moment
is quite questionable as the resulting contribution to the resistivity
due to spin fluctuations as much too small.
In fact the simulations of Ruban et al.\ \cite{RKMJ07} on the basis of a
longitudinal spin fluctuation Hamiltonian led on the case of fcc-Ni to a
strong diminishing of the average
local magnetic moment when the critical temperature is approached from below
(about 20 \% compared to $ T=0$~K).
For bcc-Fe, the change is much smaller (about 3 \%) justifying on the case the
assumption of a rigid spin moment.
Taking the extreme point of view that the spin moment vanishes completely above
the critical temperature or the paramagnetic  (PM) regime only
thermal displacements have to be considered as a source for a finite resistivity.
Corresponding results are shown in
Fig.\ \ref{FIG-RHO-EXPT} (bottom) together with corresponding
experimental data. The very good agreement between both obviously suggests that
remaining spin fluctuations above the critical temperature are of minor importance
for the resistivity of fcc-Ni.

\subsection{Temperature dependent Gilbert damping parameter}

Fig.\ \ref{FIG-GIL-THEORY} shows results for Gilbert damping parameter $\alpha$ of bcc-Fe
obtained using different models for the spin fluctuations.
All curves show the typical conductivity-like behavior
for low temperatures and the resistivity-like
behavior at high temperatures
reflecting the change from dominating
intra- to inter-band transitions.\cite{GIS07}
\begin{figure}[hbt]
 \begin{center}
   \includegraphics[angle=0,width=0.8\linewidth,clip]{Fig_8_CMP_Gilbert_fluct_Fe_lmax3.eps}
  \caption{\label{FIG-GIL-THEORY} Temperature dependent Gilbert damping
    $\alpha(T)$ for bcc-Fe, obtained by accounted for thermal vibrations
    and spin fluctuations accounting for spin fluctuations using
scheme MC (squares),
  DLM(MC) (triangles up),
 cone(MC) (triangles down)
and  a MF fit to the experimental
    temperature reduced magnetic moment (circles).
}
 \end{center}
\end{figure}
The curve denoted expt is based on a
spin configuration toted to the experimental
$ M_{ \rm expt}(T) $ data.
Using the conical model to fit $ M_{ \rm expt}(T) $
as basis for the calculation of $ \alpha (T)$
leads obviously to a rather good agreement with
$ \alpha_{ M(\rm expt) }(T)$.
Having instead a DLM-like representation of
$ M_{ \rm expt}(T) $, on the other hand,
transverse spin components are suppressed and noteworthy
deviations from $ \alpha_{ M(\rm expt) }(T)$
are found for the low temperature regime.
Nevertheless, the deviations
are less pronounced than in the case of the longitudinal
resistivity (see Fig.\ \ref{FIG-RHO-DISP-FLUCT} (bottom)), where corresponding
results are shown based on $ M_{ \rm MC}(T) $
as a reference.
Obviously, the damping parameter $ \alpha $ seems to be less sensitive to the specific spin fluctuation model
used than the resistivity.
Finally, using the spin configuration deduced from Monte Carlo simulations, i.e.\
based on $ M_{ \rm MC}(T) $
quite strong deviations for the resulting
 $ \alpha_{ M(\rm MC) }(T)$
from  $ \alpha_{ M(\rm expt) }(T)$
are found.
As for the resistivity (see Fig.\ \ref{FIG-RHO-DISP-FLUCT} (bottom))
this seems to reflect the too fast drop of the
reduced magnetic moment $ M_{ \rm MC}(T) $
 with temperature
in the low temperature regime compared
with temperature
(see Fig.\ \ref{FIG-MAG-T}).
As found before \cite{EMKK11}
accounting only for thermal vibrations
  $ \alpha_{ }(T)$ (Fig.\ \ref{FIG-RHO-DISP-FLUCT} (bottom))
is found comparable to the case when only
thermal span fluctuations are allowed.
Combing both thermal effects does not lead
to a curve that is just the sum of the two
 $ \alpha_{ }(T)$ curves.
As found for
 the conductivity (Fig.\ \ref{FIG-RHO-DISP-FLUCT} (top))
obviously the two thermal
effects are not simply additive.
As Fig.\ \ref{FIG-GIL-EXPT} (top) shows, the resulting damping parameter
 $ \alpha_{ }(T)$
for bcc-Fe
 that accounts for thermal vibrations as well as
spin fluctuations is found in reasonable good
agreement with experimental data.\cite{EMKK11}
\begin{figure}[thbt]
 \begin{center}
   \includegraphics[angle=0,width=0.8\linewidth,clip]{Fig_9a_CMP_Gilbert_fluct+vib_Fe_lmax3.eps}
   \includegraphics[angle=0,width=0.8\linewidth,clip]{Fig_9b_Ni_Gilbert_alfa_vs_T_vib_fluct.eps}
  \caption{\label{FIG-GIL-EXPT}
Top:
Temperature dependent Gilbert damping
    $\alpha(T)$ for bcc-Fe, obtained by accounted for thermal vibrations
    and spin fluctuations
 accounting for
lattice vibrations only (circles) and
    lattice vibrations and spin fluctuations based on mean-field fit to
    the experimental temperature reduced magnetic moment $M_{\rm expt} $
    (diamonds)
compared with experimental data
 (dashed and full lines).\cite{BL74,HF66}
Bottom: corresponding results
for fcc-Ni.
Experimental data have been taken from Ref.\ \onlinecite{BL74}.
}
 \end{center}
\end{figure}

 Fig.\ \ref{FIG-GIL-EXPT}
shows also corresponding results for the
Gilbert damping of fcc-Ni as a function of temperature.
Accounting only for thermal spin fluctuations on the basis of the experimental
 $M(T)$-curve leads in this case to   completely unrealistic results
while accounting only for thermal displacements leads to results already
in rather good agreement with experiment. Taking finally both sources of disorder
into account  again no simple additive behavior
is found but the results are nearly unchanged
compared to those based on the thermal displacements alone.
This implies that results for the Gilbert damping parameter of fcc-Ni
hardly depend on the specific spin configuration model used but are much
more governed by thermal displacements.

\section{Summary}
Various schemes
based on the alloy analogy that allow to include
thermal effects
when calculating
response properties relevant in spintronics
have been presented and discussed.
Technical details
of an implementation
within the framework of the spin-polarized
relativistic KKR-CPA band structure method
have been outlined that allow to
 deal with thermal
vibrations as well as spin fluctuations.
Various models to represent spin fluctuations
have been compared with each other
concerning corresponding results for the
temperature dependence of the reduced magnetic moment
$M (T)$ as well as response quantities.
It was found that response quantities
are much more sensitive to the spin fluctuation
model as the reduced magnetic moment
$M (T)$.
Furthermore, it was found that the influence
of thermal
vibrations and spin fluctuations
is not additive when calculating electrical conductivity
or the Gilbert damping parameter $\alpha $.
Using experimental data for the reduced magnetic moment
$M (T)$ to set up
 realistic temperature dependent
spin configurations satisfying agreement for
the electrical conductivity
as well as the Gilbert damping parameter
could be obtained for elemental ferromagnets
bcc-Fe  and fcc-Ni.

\begin{acknowledgments}
  This    work    was   supported    financially    by   the    Deutsche
  Forschungsgemeinschaft (DFG)
within the projects EB154/20-1, EB154/21-1 and  EB154/23-1
as well as
the
priority program SPP 1538 (Spin Caloric Transport)
 and the
SFB 689 (Spinph\"anomene in reduzierten Dimensionen).
Helpful discussions with Josef Kudrnovsk\'y and Ilja Turek
are gratefully acknowledged.
\end{acknowledgments}

redpajama_set_name: RedPajamaArXiv