## 3 Model Definition
A local density of the chosen type ions can be written, in analogy to semiclassical Pines-Bohm random phase approximation (RPA) approach to electrons in metal [2], in the following form,  
\[\rho({\bf r},t)=<\Psi_{ion}(t)|\sum_{j}\delta({\bf r}-{\bf r}_{j})|\Psi_{ion}(t)>, \tag{5}\]  
where \({\bf r}_{j}\) denotes coordinate of \(j\)-th ion. The Fourier picture of the above density has the form,  
\[\tilde{\rho}({\bf k},t)=\int\rho({\bf r},t)e^{-i{\bf k}\cdot{\bf r}}d^{3}r=< \Psi_{ion}(t)|\hat{\rho}({\bf k})|\Psi_{ion}(t)>, \tag{6}\]  
where the 'operator' \(\hat{\rho}({\bf k})=\sum\limits_{j}e^{-i{\bf k}\cdot{\bf r}_{j}}\). Using the above notation one can rewrite \(\hat{H}_{ion}\) in the following form, in an analogy to metallic plasmon description [1,2]:  
\[ \begin{array}{rl}\hat{H}_{ion}=&\sum\limits_{j=1}^{N}\left[-\frac{\hbar^{2} \
abla_{j}^{2}}{2m}\right]-\frac{q^{{}^{\prime}2}}{(2\pi)^{3}}\int d^{3}k\bar{ n}({\bf k})\frac{2\pi}{k^{2}}\left(\hat{\rho}^{+}({\bf k})+\hat{\rho}({\bf k}) \right)\\ &+\frac{q^{{}^{\prime}2}}{(2\pi)^{3}}\int d^{3}k\frac{2\pi}{k^{2}}\left[\hat{ \rho}^{+}({\bf k})\hat{\rho}({\bf k})-N\right],\end{array}  \tag{7}\]  
where: \(\bar{n}({\bf k})=\int d^{3}rn({\bf r})e^{-i{\bf k}\cdot{\bf r}}\) (we have taken into account that \(\frac{4\pi}{k^{2}}=\int d^{3}r\frac{1}{r}e^{-i{\bf k}\cdot{\bf r}}\), \(q^{{}^{\prime}2}=\frac{q^{2}}{\varepsilon}\)). Utilizing this form of the effective ion Hamiltonian one can write out the dynamic equation in Heisenberg representation for ion density fluctuations,  
\[\frac{d^{2}\hat{\rho}({\bf k})}{dt^{2}}=\frac{1}{(\bar{n})^{2}}\left[\left[ \hat{\rho}({\bf k}),\hat{H}_{ion}\right],\hat{H}_{ion}\right], \tag{8}\]  
in the following form,  
\[ \begin{array}{rl}\frac{d^{2}\hat{\rho}({\bf k})}{dt^{2}}=-\sum\limits_{j}e^{ -i{\bf k}\cdot{\bf r}_{j}}\left\{-\frac{\hbar^{2}}{m^{2}}\left({\bf k}\cdot \
abla_{j}\right)^{2}+\frac{\hbar^{2}k^{2}}{m^{2}}i{\bf k}\cdot\
abla_{j}+\frac {\hbar^{2}k^{4}}{4m^{2}}\right\}\\ -\frac{4\pi q^{{}^{\prime}2}}{m(2\pi)^{3}}\int d^{3}p\bar{n}({\bf p})\frac{{ \bf k}\cdot{\bf p}}{p^{2}}\hat{\rho}({\bf k}-{\bf p})-\frac{4\pi q^{{}^{\prime} 2}}{m(2\pi)^{3}}\int d^{3}p\hat{ \rho}({\bf k}-{\bf p})\frac{{\bf k}\cdot{\bf p}}{p^{2}}\hat{\rho}({\bf p}). \end{array}  \tag{9}\]  
If one takes into account that \(\hat{\rho}({\bf k}-{\bf p})\hat{\rho}({\bf p})=\delta\hat{\rho}({\bf k}-{\bf p })\delta\hat{\rho}({\bf p})+\bar{n}({\bf k}-{\bf p})\delta\hat{\rho}({\bf p})+ \delta\hat{\rho}({\bf k}-{\bf p})\bar{n}({\bf p})+\bar{n}({\bf k}-{\bf p})\bar {n}({\bf p})\) and \(\bar{n}({\bf p})\hat{\rho}({\bf k}-{\bf p})=\bar{n}({\bf p})\delta\hat{\rho}( {\bf k}-{\bf p})+\bar{n}({\bf p})\bar{n}({\bf k}-{\bf p})\), where \(\delta\hat{\rho}({\bf k})=\hat{\rho}({\bf k})-\bar{n}({\bf k})\) describes the 'operator' of local ion density fluctuation above the uniform distribution, one can rewrite Equation (9) as follows,  
\[ \begin{array}{rl}\frac{d^{2}\delta\hat{\rho}({\bf k})}{dt^{2}}=-\sum\limits_{j }e^{-i{\bf k}\cdot{\bf r}_{j}}\left\{-\frac{\hbar^{2}}{m^{2}}\left({\bf k} \cdot\
abla_{j}\right)^{2}+\frac{\hbar^{2}k^{2}}{m^{2}}i{\bf k}\cdot\
abla_{j}+ \frac{\hbar^{2}k^{4}}{4m^{2}}\right\}\\ -\frac{4\pi q^{{}^{\prime}2}}{m(2\pi)^{3}}\int d^{3}p\bar{n}({\bf k}-{\bf p}) \frac{{\bf k}\cdot{\bf p}}{p^{2}}\delta\hat{\rho}({\bf p})-\frac{4\pi q^{{}^{ \prime}2}}{m(2\pi)^{3}}\int d^{3}p\delta\hat{\rho}({\bf k}-{\bf p})\frac{{\bf k }\cdot{\bf p}}{p^{2}}\delta\hat{\rho}({\bf p}).\end{array}  \tag{10}\]  
Hence, for the ion density fluctuation, \(\delta\hat{\rho}({\bf k},t)=<\Psi_{ion}|\delta\hat{\rho}({\bf k},t)|\Psi_{ion}>= \hat{\rho}({\bf k},t)-\bar{n}({\bf k})\), one can find,  
\[