in Sec. \[Sec:electron-correlation-eff\]. Finally, we estimate theoretical accuracy and consider its implications on the extraction of octupole and hexadecapole moments in Sec. \[sec:Discussion\]. Unless specified otherwise, atomic units are used throughout.

Review of the theory of hyperfine structure {#Sec:General-theory}
===========================================

The hyperfine interaction can be decomposed into the magnetic dipole (M1), electric quadrupole (E2), magnetic octupole (M3), electric hexadecapole (E4), and higher rank contributions. We start by expressing the hyperfine Hamiltonian in irreducible tensor form [@sch55; @Johnson2007] $$H_\mathrm{HFI}=\sum\limits_{k,\mu}(-1)^{\mu}T_{k,\mu}^eT_{k,-\mu}^n \,,\ \label{Eq:HF-general}$$ where rank-$k$ tensors $T_{k,\mu}^e$ act in the electron space, and $T^n_{k,-\mu}$ — in the nuclear space. The many-electron operators are $T_{k,\mu}^e=\sum\limits_{i}t_{k,\mu}^e(i)$, where the summation is over all the atomic electrons. The single-electron operators $t_{k,\mu}^e(i)$ can be divided into two groups [@Johnson2007], $$t_{k,\mu}^e(i)= \begin{cases}
      -\frac{1}{r^{k+1}}C_{k,\mu}(\hat{r}), & \text{electric (even $k$)}\,, \\
      -\frac{i}{r^{k+1}}\sqrt{\frac{k+1}{k}} \bm{\alpha} \cdot{\bm{C}^{(0)}_{k,\mu}}(\hat{r}), & \text{magnetic (odd $k$)}. 
      \end{cases}\label{Eq:single-electron-operator}$$ Here, $\bm{\alpha}$ is the Dirac matrix, $r$ is the radial coordinate, $C_{k,\mu}$ are normalized spherical harmonics, and $\bm{C}_{k,\mu}^{(0)}$ are normalized vector spherical harmonics.

The first-order energy correction due to hyperfine interaction, Eq. (\[Eq:HF-general\]), in the basis of coupled nuclear and atomic states is [@Beloy2008a] $$\begin{split}
 W^{(1)}_F = \langle \gamma{IJFM_F}|H_{\mathrm{HFI}}|\gamma IJFM_F\rangle = (-1)^{I+J+F} \\
    \times\sum\limits_{k}\begin{Bmatrix}
    F & J & I \\
    k & I 