## 1 Introduction
### Nonlinear Effects Under Ion-Relaxation Polarization Conditions  
A phenomenological model of the diffusion transfer of ions in ion dielectrics (in the case of HBCs, protons) in an electric field, linearized in the generalized kinetic Equation (28) via the parameter k (with k = 1 in (40a)), was developed on the basis of a system of nonlinear equations of the Fokker-Planck and Poisson types [1, 122, 133, 135]:  
\[\frac{\partial\rho}{\partial\tau}=\frac{\partial^{2}\rho}{\partial\xi^{2}}- \theta\rho-\gamma\frac{\partial}{\partial\xi}(\rho\mathsf{z}), \tag{57}\]  
\[\frac{\partial\mathsf{z}}{\partial\xi}=\phi_{\mathrm{I}}\rho, \tag{58}\]  
and their initial and boundary conditions are  
\[\rho(\xi,0)=0, \tag{59}\]  
\[\frac{\partial\rho}{\partial\xi}\bigg{|}_{\xi=\{0;\frac{\mathrm{d}}{2}\}}= \gamma(n_{0}+\rho)\,\mathsf{z}\big{|}_{\xi=\{0;\frac{\mathrm{d}}{2}\}}, \tag{60}\]  
\[\int_{0}^{\frac{\mathrm{d}}{2}}z(\xi,\tau)\mathrm{d}\xi=\frac{\mathrm{d}}{a} \exp\biggl{(}\frac{\mathrm{i}\omega}{W^{(0)}}\tau\biggr{)}. \tag{61}\]  
In (57)-(61), the following designations are adopted: \(\rho(\xi,\tau)=n(\mathsf{x};\mathsf{t})-n_{0}\) denotes the concentration of ions excessive over their equilibrium concentration \(n_{0}\), \(\varepsilon_{\infty}\) is the high-frequency dielectric constant of the crystal, \(\mathsf{z}(\xi,\tau)=\frac{\mathsf{E}(\mathsf{x}\tau)}{\mathsf{E}_{0}}\), \(\xi=\frac{\mathsf{x}}{a}\), \(\tau=W^{(0)}\mathsf{t}\), \(\phi_{\mathrm{I}}=\frac{\omega\mathsf{z}}{\varepsilon_{0}\varepsilon_{0} \varepsilon_{0}}\), and \(\theta=\phi_{\mathrm{I}}\gamma n_{0}\), \(\gamma=\frac{\varepsilon_{0}W^{(0)}}{W^{(0)}}\). A solution to the system of Equations (3) and (4) deriving from perturbation theory methods was developed using the following power series [1, 125, 135]:  
\[\rho(\xi;\tau)=\sum_{k=1}^{+\infty}\gamma^{k}\rho_{\mathrm{k}}(\xi;\tau), \mathsf{z}(\xi;\tau)=\sum_{k=0}^{+\infty}\gamma^{k}\mathsf{z}_{\mathrm{k}}( \xi;\tau). \tag{62}\]  
As shown in the Appendix A, based on expressions (A9a)-(A9j), the relaxation modes \(\rho_{\mathrm{k}}^{(\omega)}(\xi,\tau)\) of the components of the volumetric charge density \(\rho(\xi,\tau)\) calculated in the k-th approximation of perturbation theory at the frequencies \((\omega;2\omega)\) of the alternating field (expressions (A10a) and (A10b)) were determined via mathematical induction. An attempt to apply a method similar to the calculation of \(\rho(\xi;\tau)\) in the k approximation at frequency \(3\omega\) yields such cumbersome expressions, such as \(\rho_{4}^{(3\omega)}(\xi,\tau)\), \(\rho_{5}^{(3\omega)}\), etc., that the derivation of the recurrence expression \(\rho_{\mathrm{k}}^{(3\omega)}(\xi,\tau)\) requires a fundamentally different and more general approach. The simplest is the expression \(\rho_{3}^{(3\omega)}(\xi,\tau)\) formulated in (A9f). Then, by substituting series (62) into the system (57)-(61), we obtain  
\[\frac{\partial\rho_{\mathrm{k}}}{\partial\tau}=\frac{\partial^{2}\rho_{ \mathrm{k}}}{\partial\xi^{2}}-\frac{\partial}{\partial\xi}\Bigl{(}\mathsf{z} _{0}\rho_{\mathrm{k}-1}+\sum_{m=1}^{\mathrm{k}-2}\mathsf{z}_{m}\rho_{\mathrm{ k}-m-1}\Bigr{)}-\theta\rho_{\mathrm{k}}, \tag{63}\]  
\[\frac{\partial\mathsf{z}_{\mathrm{k}}}{\partial\xi}=\phi\rho_{\mathrm{k}}, \tag{64}\]  
\[\rho_{\mathrm{k}}(\xi;0)=0, \tag{65}\]  
\[\frac{\partial\rho_{\mathrm{k}}}{\partial\xi}\bigg{|}_{\xi=(0;\frac{\mathrm{ d}}{2})}=\left.\Bigl{[}n_{0}\mathsf{z}_{\mathrm{k}-1}+\mathsf{z}_{0}\rho_{ \mathrm{k}-1}+\sum_{m=1}^{\mathrm{k}-2}\mathsf{z}_{m}\rho_{\mathrm{k}-m-1} \Bigr{]}\right|_{\xi=(0;\frac{\mathrm{d}}{2})}. \tag{66}\]  
Here, \(\rho_{0}(\xi;\tau)=0\), and \(\mathsf{z}_{0}(\xi;\tau)=\exp\Bigl{(}\frac{\mathrm{i}\omega}{W^{(0)}}\tau \Bigr{)}\)[1, 3, 5].