fitting is good enough to reproduce results from MC simulations with parameters shown in and Tab. \[tab:parameters-BZT-static\]. The closeness of $E_{b}$ and $T_{m}$ indicates the average depth of potential wells plays a dominant role in determining $T_{m}$. Close examination of Fig. \[fig:Fitting-susceptibility-of\] also indicates the fitting at the lowest temperature ($\lesssim25$ K) is not as good as the rest. To address this issue, we tried adding a Gaussian distribution to $E_{b}$ and remedied the minor problem. However, the resulting equation is quite complicated and deviates from our original goal of proposing simple analytical formula to fit susceptibility. Therefore this additional step is not adopted here.

            $\chi_{1}$   $\theta$ (K)   $\chi_{2}$   $E_{b}$ (K)
  -------- ------------ -------------- ------------ -------------
   Values     741.6         220.5          64.7         159.1

  : Fitting parameters for the Ba(Zr$_{0.5}$,Ti$_{0.5}$)O$_{3}$ static susceptibility. \[tab:parameters-BZT-static\]

![Maxwell-Boltzmann distribution [\[]{}Panel (a)[\]]{} and the Langevin function, $\mathcal{L}$, [\[]{}Panel (b)[\]]{} versus temperature. Parameters from Tab. \[tab:parameters-BZT-static\] are used in plot each function. \[fig:components\]](Fig_0){width="6cm"}

In order to have a good understanding of BZT’s susceptibility, we show each component of Eq. (\[eq:BZT-susceptibility\]) in Fig. \[fig:components\]. Figure \[fig:components\] (a) shows that $P_{1}\left(E_{b},T\right)$ and $P_{2}\left(E_{b},T\right)$ have opposite trends as temperature increases. The number of dipoles that can overcome the potential confinement ($P_{1}$) steadily 