## 1 Introduction
### Complex Dielectric Permittivity
By assuming that \(\lambda=0\) when \(\Omega_{1}=\{\omega\}\) and \(\hat{\alpha}_{0}^{(\Omega_{1}=\{\omega\})}=\hat{\alpha}^{(\Omega_{1}=\{\omega \})}=\frac{\mathrm{P}^{(\omega)}(\mathrm{t})}{\mathrm{E}_{\mathrm{pol}}( \mathrm{t})}\), owing to (A14a), we have  
\[\hat{\alpha}_{0}^{(\Omega_{1}=\{\omega\})}=\hat{\alpha}^{(\Omega_{1}=\{\omega \})}=\frac{\varepsilon_{0}\varepsilon_{\infty}\Big{(}\Gamma_{1}^{(\omega)}- \mathrm{i}\Gamma_{2}^{(\omega)}\Big{)}}{1-\Gamma_{1}^{(\omega)}+\mathrm{i} \Gamma_{2}^{(\omega)}}, \tag{99a}\] \[\hat{\varepsilon}_{0}^{(\Omega_{1}=\{\omega\})}=\frac{\hat{\alpha}_{0}^{(\Omega _{1}=\{\omega\})}}{\varepsilon_{0}}=\frac{\mathrm{P}^{(\omega)}(\mathrm{t})}{ \varepsilon_{0}\mathrm{E}_{\mathrm{pol}}(\mathrm{t})}=\frac{\varepsilon_{ \infty}\Gamma^{(\omega)}}{1-\Gamma^{(\omega)}}. \tag{99b}\]  
In (99a) and (99b), due to the great complexity of this research, corresponding to \(\Gamma^{(\omega)}=\Gamma_{1}^{(\omega)}-\mathrm{i}\Gamma_{2}^{(\omega)}\), we use the infinite sums of ranks, and the material and imaginary components of the parameter \(\Gamma^{(\omega)}\) are considered in the kinetic theory of polarization to be quasiclassical relaxation parameters that are multiples of the main frequency \((\omega)\). \[\mathrm{P}^{((2\lambda+1)\omega)}(\mathrm{t})=\frac{\mathrm{q}}{\mathrm{d}} \mathrm{\sum}_{\lambda=0}^{\infty}\int_{0}^{\mathrm{d}}\mathrm{x}\rho^{((2 \lambda+1)\omega)}(\xi,\tau)\mathrm{d}\mathrm{x},\]0  
In (100), \(T_{n}=\frac{\tau_{n}}{W^{(0)}}=\Big{(}\frac{1}{T_{n,\mathrm{D}}}+\frac{1}{T_{ \mathrm{M}}}\Big{)}^{-1}\) denotes the relaxation time for the n-th relaxation mode, i.e., \(T_{n}=\frac{\tau_{n,\mathrm{D}}T_{\mathrm{M}}}{T_{n\mathrm{D}}+T_{\mathrm{M}}}\), where \(T_{n,\mathrm{D}}=\frac{T_{\mathrm{D}}}{\mathrm{n}^{2}}\) is the diffusion relaxation time for the n-th relaxation mode and \(T_{\mathrm{D}}=\frac{\mathrm{d}^{2}}{\pi^{2}\mathrm{D}_{\mathrm{diff}}^{(0)}}\) is that for the 1 st relaxation mode. \(T_{\mathrm{M}}=\frac{\varepsilon_{0}\varepsilon_{\infty}}{\mathrm{u}_{\mathrm{ mD}}^{(1)}\cdot\mathrm{q}\;\mathrm{n}_{0}}\) denotes the Maxwell relaxation time. As noted in Section 2.5, the coefficients \(\mathrm{D}_{diff}^{(0)}=a^{2}\cdot W^{(0)}\) and \(\mathrm{u}_{\mathrm{mD}}^{(1)}=\frac{\mathrm{q}a^{2}W^{(1)}}{k_{\mathrm{B}}T}\) are computed at \(l=0\) in the nonlinear coefficients \(\mathrm{D}_{diff}^{(2l)}=a^{2}\cdot W^{(2l)}\) and \(\mathrm{u}_{\mathrm{mD}}^{(2l+1)}=\frac{\mathrm{q}a^{2}W^{(2l+1)}}{k_{\mathrm{B}} T}\) from (19). In (100), \(\tau_{n}\) is the dimensionless relaxation time for n relaxation modes. Furthermore, moving on to the following order of calculations, assuming \(\lambda=1\) when \(\Omega_{3}=\{\omega;2\omega;3\omega\}\) and \(\hat{\alpha}^{(\Omega_{3}=\{\omega;2\omega;3\omega\})}=\hat{\alpha}_{2}^{( \Omega_{3}=\{\omega;2\omega;3\omega\})}\cdot\mathrm{E}_{\mathrm{pol}}^{2 \lambda}(\mathrm{t})=\frac{\mathrm{P}^{(3\omega)}(\mathrm{t})}{\mathrm{E}_{ \mathrm{pol}}(\mathrm{t})}\) from (97), owing to (92), we have  
\[\mathrm{P}^{((2\lambda+1)\omega)}(\mathrm{t})=\frac{\mathrm{q}}{\mathrm{d}} \mathrm{\sum}_{\lambda=0}^{\infty}\int_{0}^{\mathrm{d}}\mathrm{x}\rho^{((2 \lambda+1)\omega)}(\xi,\tau)\mathrm{d}\mathrm{x},\]1  
Accordingly,  
\[\mathrm{P}^{((2\lambda+1)\omega)}(\mathrm{t})=\frac{\mathrm{q}}{\mathrm{d}} \mathrm{\sum}_{\lambda=0}^{\infty}\int_{0}^{\mathrm{d}}\mathrm{x}\rho^{((2 \lambda+1)\omega)}(\xi,\tau)\mathrm{d}\mathrm{x},\]2Here, \(\Upsilon=\frac{\varepsilon_{0}^{2}\varepsilon_{\infty}^{2}\pi^{4}}{4\Omega^{2} \eta_{0}^{2}\Omega^{2}}\).