How is the two-level atom dipole force modified for a real, multi-level atom? I was trying to compute the dipole force exerted by a laser beam at $1064$ nm on a $^6{\mathrm{Li}}$ atom (the most important transition being around $671$ nm, the laser is greatly red-detuned), when I came upon a conceptual difficulty.
I started from the usual formula for a two-level atom (see Grimm 1999 for instance):
$$U_{\mathrm{dip}}(r) = - \frac{3 \pi c^2}{2 \omega_0^3} \left(\frac{\Gamma}{\omega_0 - \omega} + \frac{\Gamma}{\omega_0 + \omega}\right) I(r),$$
where $\Gamma$ is the linewidth of the considered transition, $\omega_0$ its frequency and $I(r)$ is the local light intensity. Now I have nothing wrong with this formula itself, and I roughly know how to derivate it. However, I was wondering how does this extend to real, multi-level atoms?
I know for instance that this relies partly on the calculation of $\Gamma$ for a two-level atom:
$$\Gamma = \frac{\omega_0^3}{3 \pi \varepsilon_0 \hbar c^3} |d|^2,$$ with $d = \left\langle e \right| q \hat{r} \left| g \right\rangle$ the dipole matrix element of the transition. But I am pretty sure that this only holds for a two-level system, otherwise you have to sum over the different possible final states.
For a real atom, like $^6{\mathrm{Li}}$, I would guess that you would have to sum something over all possible excited states, possibly using a formula that lets explicitly appear $d$ instead of $\Gamma$. If you consider an atom in its ground state $2^2 S_{1/2}$, it is mostly connected to two excited manifolds, which are $2^2 P_{1/2}$ (called the $D_1$ line), and $2^2 P_{3/2}$ ($D_2$ line) (see Gehm 2003). This is already complicated, but these two manifolds can be decomposed respectively into $2$ and $3$ additional manifolds corresponding to different $F$ values. Do you have to sum over all of these levels, even though the beam is clearly commpletely off-resonant so I would guess it "doesn't care" about the exact hyperfine structure of $^6{\mathrm{Li}}$.
To even add to my confusion, my advisor just told me: "Oh, but look! $D_1$ and $D_2$ almost have the same $\Gamma \approx 2\pi \times 6$ MHz, so just plug $\Gamma = 2\pi \times 6$ MHz in the above formula." Aren't you supposed to add them instead (and end up with roughly a factor $2$)?
I am as interested in getting an "efficient" way for us, experimentalists, to compute the dipole potential for a given atom, knowing all the transitions linewidth/dipole elements, as I am in getting a better undestanding of the underlying calculations. Can you "forget" about the fine and hyperfine structure when considering very detuned driving? Are there any resources specifically adressing these type of issues (extending two-level atoms results to real atoms)?
Thank you very much for your help.
 A: This has been done by many people. It is called the ac-Stark shift, and maybe you are able find programs in the internet, which are already  for Lithium. Check out the phd theses in the ultra-cold atom community.
I have done it for an atom which possesses only a $(LS)J$ coupling. In your case, you will have to adapt the methods, which will be painful -- I know what I'm talking about, because I once had do exactly the same calculate for Lithium as well. Unfortunately, I don't have these calculations anymore. 
Using the polarisation 
$\vec \epsilon_{q}$, the ac-Stark shift
of the state $|\alpha\rangle$ is
given by
$$
    \Delta E_{\alpha}^{(q)} = %
    \frac{E_{0}^2}{4 \pi \hbar c} \sum_{\alpha^{\prime}} \;%
    sign_{\alpha^{\prime}, \alpha} \;%
    % \left| \; %
    |\mu_{\alpha^{\prime}, q, \alpha}|^2 %
    \frac{\lambda_{\alpha^{\prime}, \alpha} \, %
      \lambda_{\textrm{laser}}^2}{%
      \lambda_{\textrm{laser}}^2%
      - \lambda_{\alpha^{\prime}, \alpha}^2} %
    % \; \right|%
$$
with the atomic transition wavelength $\lambda_{\alpha^{\prime}, \alpha}$ and the (electric dipole) transition matrix element
$$
\mu_{\alpha^{\prime}, \vec q,  \vec \alpha} =%
    e \, \langle\alpha^{\prime}| %
    r \cdot \epsilon_{q}|\alpha\rangle
%%%
= e\, 
\langle \gamma^{\prime} n^{\prime} (S L^{\prime}) J^{\prime} 
m_{J^\prime} |
\vec r \cdot \vec \epsilon_{q} 
 |\gamma n (S L) J m_{\textrm{J}} \rangle
$$ 
Hint: The sign of the ac-Stark shift should be put explicitly.
  $sign_{\alpha^{\prime}, \alpha}=-1$ for the energetically
  lower state, while it is $+1$ for the upper state. 
Now, what I did is to use the NIST database to obtain all atomic levels and all coupling strengths. You will have to download this and bring it into a format, so that you can use it in your program. This took a lot of work, however, it was almost 10 years ago. Maybe things have changed. 
Be aware that NIST defines 
$$
|\mu_{\alpha^{\prime}, q, \alpha}|^2 
=
\frac{3 \, \epsilon_{0}\, h \,
      \lambda_{\alpha^\prime, \alpha}^3}{16 \pi^3} \, %
    A^{(\textrm{NIST})}_{\vec \alpha^\prime, \alpha} %
    \times %
    \langle{\gamma n J \, 1\, m_{\textrm{J}} \,
      q|}{\gamma^{\prime} n^{\prime} (J 1)J^{\prime} \,
      m_{\textrm{J}^{\prime}}}\rangle^{2}
$$
where the last term is the square of the Clebsch-Gordan coefficient. Since Lithium 
possesses a nuclear spin $I_{\textrm{spin}}$ the formula
is given by  (please cross-check)
\begin{align}
  \label{eq:ACStarkShift}
  E_{j} = \frac{3\pi}{2}\, c^2  \, I \sum_{k} \;%
  &\frac{A_{k,j}}{\omega_{k,j}^3} %
  \left(%
    \frac{1}{\omega_{k,j}+\omega_{\textrm{laser}}} + % NO MINUS !!!
    \frac{1}{\omega_{k,j}-\omega_{\textrm{laser}}}%
  \right) \\%
  &\times %
  (2\,J_k +1)   \, %
  (2\,F_j+1) \, (2\,F_k+1) %
  \nonumber \\%
  &\times %
  \left|%
    \begin{pmatrix} %
      J_k & 1 & J_j \\ %
      -M_k & q & M_j %
    \end{pmatrix} %
    \left\{ %
    \begin{matrix} %
      J_k & I_{\textrm{spin}} & F_k \\ %
      F_j & 1 & J_j %
    \end{matrix} %
    \right\} %
  \right|^2 %
  \nonumber %
\end{align}
  Notation: The  indices $j$ and $k$ are for the energetically lower
  and upper state, %are understood to be multi-indices, 
  $q$ specifies the polarisation of the
  light ($q=0, \pm1$ for $\pi$, $\sigma^{(\pm)}$ light), 
  and the round $(\ldots)$ and curly $\{\ldots\}$ brackets denote the
  Wigner-3j and Wigner-6j symbols, respectively. 
These things always look straight forward, however, once you start to implement it, you quickly get confused. So my advice would be that you take a piece of paper and write done all equations to a level, which seams to be ridicules. Follow each step through, and verify each equation. Almost every paper and phd thesis which included this topic, had several misprints or short-comings.  Good luck, and don't get too frustrated. 
