buckling. The $R_B$ ($=2.732$Å) and $R_C$ ($=2.676$Å) are distances from the origin to B and C sites, respectively. The $d_B$ and $d_C$ represent the absolute values of shifts along the $c$ axis from the $z=\pm c/2$ planes, respectively. Angles $\theta_{B(C)}^{(j)}$ and $\theta_{B(C)}^{(j)}$ in eqs. (\[eq.defA22\]) and (\[eq.defAxy\]) are given by $$\begin{aligned}
 \theta_B^{(1)} &= 65.6^{\circ}, \quad \theta_C^{(1)} = 19.1^{\circ},\nonumber\\
 \theta_B^{(2)} &= 180^{\circ}-65.6^{\circ}, 
   \quad \theta_C^{(2)} = 180^{\circ}-19.1^{\circ},\nonumber\\
 \theta_B^{(3)} &= 90^{\circ}+\theta_B^{(1)}, 
   \quad \theta_C^{(3)} = 90^{\circ}+\theta_C^{(1)},\nonumber\\
 \theta_B^{(4)} &= 90^{\circ}+\theta_B^{(2)}, 
   \quad \theta_C^{(4)} = 90^{\circ}+\theta_C^{(2)}.\end{aligned}$$

Now we search for the local coordinate frames in which the third term in eq. (\[eq.Vcry\]) is eliminated. Rotating the original coordinate frame by angle $\phi_j$ around the $c$ axis for each sublattice, we have the operators transformed as $$\begin{aligned}
 Q_2^2 &= \cos(2\phi_j)\tilde Q_2^2(j) - \sin(2\phi_j)\tilde Q_{xy}(j),
\nonumber\\
 Q_{xy} &= \sin(2\phi_j)\tilde Q_2^2(j) + \cos(2\phi_j)\tilde Q_{xy}(j),
\label{eq.rotQ}\end{aligned}$$ where tilde operators $\tilde Q(j)$’s are represented with respect to the local coordinate frames. $Q_2^0$ is unchanged. Inserting eq. (\[eq.rotQ\]) into eq. (\[eq.Vcry\]), we have $$V_{crys}(j) = A_2^0 \tilde Q_2^0(j) + \tilde A_2^2(j)\tilde Q_2^2(j)
             + \tilde A_{xy}(j)\tilde Q_{xy}(j),$$ with $$\begin{aligned}
 \tilde A_{2}^{2}(j) &= \Lambda_{B} \cos [2(\theta_{B}^{(j)}-\phi_j)]
            + \Lambda_{C} \cos [ 2 (\theta_{C}^{(j)}-\phi_j)], \nonumber \\
\\
 \tilde A_{xy}(j) &= \Lambda_{B} \sin [2(\theta_{B}^{(j)}-\phi_j)]
          + \Lambda_{C} \sin [2(\theta_{C}^{(j)}-\phi_j)]. 