We derive explicit expressions for the magnetic susceptibility at high temperatures for the full Hamiltonian.
Including exchange interactions, our full Hamiltonian reads
\[\begin{split}
{\cal H}
&=
\sum_i{\cal H}_\text{cef}(i)
+\sum_{\left\langle ij\right\rangle}{\cal H}_\text{exc}(ij)
+\sum_i{\cal H}_\text{Zeeman}(i),
\\
{\cal H}_\text{cef}(i)
&=
\sum_{\ell=2\atop\Delta\ell=2}^6\sum_{m=0\atop\Delta m=3}^\ell
B_\ell^mO_\ell^m(i),
\\
{\cal H}_\text{exc}(ij)
&=
\sum_{\alpha\beta}J_i^\alpha J_{ij}^{\alpha\beta}J_j^\beta,
\\
{\cal H}_\text{Zeeman}(i)
&=
-g_j\mu_\text B\mu_0\sum_\alpha H_i^\alpha J_i^\alpha.
\end{split}\]
The dimensionless uniform magnetic susceptibility $\chi(T)$ of a crystal with volume $V$ is given by the change of the magnetisation $M$ with the magnetic field $B=\mu_0H$ with components
\[\chi_\alpha=\mu_0\frac{\partial M_\alpha}{\partial B_\alpha}
=\frac{\mu_0}V\frac{\partial \bar\mu_\alpha}{\partial B_\alpha}
=-\frac{\mu_0}V\frac{\partial^2F}{\partial B_\alpha^2},
\quad\alpha=\parallel,\perp\]
where $\bar\mu_\alpha=-\partial F/\partial B_\alpha$ is the total magnetic moment in spatial direction $\alpha$ either parallel or perpendicular to the c axis and $F=(1/\beta)\log\cal Z$ the canonical free energy of the Hamiltonian above, $1/\beta=k_\text BT$ the inverse temperature and $k_\text B$ the Boltzmann constant. For the molar susceptibility the equation has to be multiplied by $N_\text L/(\nu/V)$ where $N_\text L$ is Avogadro’s number and $\nu/V$ the volume density of $\text{Yb}^{3+}$ ions.
An expansion of $\chi_\alpha$ in powers of $\beta$, together with the orthogonality and tracelessness of the Stevens operators $O_\ell^m$ eventually yields
\[\begin{split}
k_\text B\Theta_\parallel
&=
-\frac45\left(j-\frac12\right)\left(j+\frac32\right)B_2^0
-\frac13j(j+1)J_\parallel(0),
\\
k_\text B\Theta_\perp
&=
\frac25\left(j-\frac12\right)\left(j+\frac32\right)B_2^0
-\frac13j(j+1)J_\perp(0)
\end{split}\]
with $J_\parallel(0)=6\hat J_\parallel$, $J_\perp(0)=\hat 6J_\perp$ in our case. Remarkably, due to the orthogonality and tracelessness of the Stevens operators, only the $B_2^0$ CEF parameter enters the Curie-Weiss temperatures, independent of the form and symmetry of the crystal field otherwise as long as it contains a more-than-2-fold symmetry axis.
The hat here should indicate the exchange parameters refering to the full $j=7/2$ moment. At low enough temperatures, only the pseudospin doublet enters the susceptibility. Here we obtain
\[k_\text B\Theta_\alpha=-\frac32J_\alpha,\quad\alpha=\parallel,\perp\]
with $J_\alpha=(g_\alpha/g_j)^2\hat J_\alpha$. (See this page for a definition of the g factors.)