Dy sites. As shown in Fig. \[fig.cryst\](b), four inequivalent sites arises. For sites $j$ ($=1\cdots 4$), the electrostatic potential is obtained as $$V_{\rm crys}(j) = A_2^0 Q_2^0
              + A_2^2(j) Q_2^2
              + A_{xy}(j)Q_{xy} ,
\label{eq.Vcry}$$ with $$\begin{aligned}
 Q_2^0 &= \frac{1}{2}(3z^2-r^2), \\
 Q_2^2 &= \frac{\sqrt{3}}{2}(x^2-y^2), \\
 Q_{xy}&= \sqrt{3}xy,\end{aligned}$$ where the coefficients are given by $$\begin{aligned}
A_{2}^{0} & = 
 \frac{2q_D}{a^3}
\left[ 1 -2 \frac{|q_B|}{q_D}
 \left(\frac{a}{R_{B}} \right)^3
\left\{ 1 - 3 \left( \frac{c}{2 R_{B}} \right)^2 \right\}
\right. \nonumber \\
&  \left.
         -2 \frac{|q_C|}{q_D}
  \left(\frac{a}{R_{C}} \right)^3
\left\{ 1 - 3 \left( \frac{c}{2 R_{C}} \right)^2 \right\} \right]
\nonumber \\
 &= \frac{2q_D}{a^3}\left[ 1 + 1.40\frac{|q_B|}{q_D} + 1.79\frac{|q_C|}{q_D}
  \right], \\
A_{2}^{2}(j) &= \Lambda_{B} \cos ( 2 \theta_{B}^{(j)} )
            + \Lambda_{C} \cos ( 2 \theta_{C}^{(j)} ), \label{eq.defA22} \\
A_{xy}(j) &=  \Lambda_{B} \sin ( 2 \theta_{B}^{(j)} )
          + \Lambda_{C} \sin ( 2 \theta_{C}^{(j)} ), \label{eq.defAxy}\end{aligned}$$ with $$\begin{aligned}
\Lambda_{B} & = - \frac{|q_B|}{a^3}
 \frac{60}{\sqrt{3}}  \left(\frac{r_B}{a} \right)^2
\left(\frac{a}{R_{B}}\right)^5
\left( \frac{c}{2 R_{B}} \right)^2 \frac{d_B}{\frac{c}{2}}
\nonumber \\
&= -22.3 \frac{|q_B|}{a^3} \frac{d_B}{\frac{c}{2}}, \\
\Lambda_{C} & = - \frac{|q_C|}{a^3}
 \frac{60}{\sqrt{3}}  \left(\frac{r_C}{a} \right)^2
\left(\frac{a}{R_{C}}\right)^5
\left( \frac{c}{2 R_{C}} \right)^2  \frac{d_C}{\frac{c}{2}}
\nonumber \\
&= -24.0 \frac{|q_C|}{a^3}  \frac{d_C}{\frac{c}{2}}.\end{aligned}$$ The first term in eq. (\[eq.Vcry\]) represents the crystal field without lattice distortion, while the second and third terms arise from the 