## Appendix A Derivation of Plasmon Frequencies for Finite Electrolyte System
### Surface Ionic Plasmons  
In order to determine self-frequencies for surface plasmons, one has to consider Equation (19) and the solution for it given by Equation (20). The first term in r.h.s. of Equation (19) can be rewritten to the form,  
\[ \begin{array}{l}\frac{kT}{m}\
abla(n+\delta\rho_{2})\
abla\Theta(a-r)+\frac{ k_{B}T}{m}(n+\delta\rho_{2})\Delta\Theta\\ =-\frac{k_{B}T}{m}\delta(a-r)\frac{\partial}{\partial r}(n+\delta\rho_{2})- \frac{k_{B}T}{m}(n+\delta\rho_{2})\frac{1}{r^{2}}\frac{\partial}{\partial r}( r^{2}\delta(a-r))\\ =-\frac{k_{B}T}{m}\frac{1}{r^{2}}\frac{\partial}{\partial r}\left[(n+\delta\rho _{2})r^{2}\delta(a-r)\right],\end{array}  \tag{11}\]  
where we have used the formulae \(\
abla\Theta(a-r)=-\frac{r}{r}\delta(a-r)\), \(\frac{r}{r}\
abla=\frac{\partial}{\partial r}\). The next term in r.h.s. of Equation (19) can be transformed into,  
\[ \begin{array}{l}-\frac{k_{B}T}{m}\delta(a-r)\frac{r}{r}\
abla\delta\rho_{2} -\frac{\omega_{F}^{2}}{4\pi}\delta(a-r)\frac{r}{r}\
abla\int\frac{d^{2}r_{1} \delta\rho(\mathbf{r}_{1})}{|\mathbf{r}-\mathbf{r}_{1}|}\\ =-\frac{k_{B}T}{m}\delta(a-r)\frac{\partial}{\partial r}\delta\rho_{2}-\frac{ \omega_{F}^{2}}{4\pi}\delta(a-r)\frac{\partial}{\partial r}\int\frac{d^{2}r_{1 }\delta\rho(\mathbf{r}_{1})}{|\mathbf{r}-\mathbf{r}_{1}|}.\end{array}  \tag{12}\]  
Equation (19) attains thus the form,  
\[ \begin{array}{l}\frac{\partial^{2}\rho_{2}}{\partial t^{2}}=-\frac{k_{B}T}{ m}\frac{1}{r^{2}}\frac{\partial}{\partial r}\left[(n+\delta\rho_{2})r^{2} \delta(a-r)\right]\\ -\frac{k_{B}T}{m}\delta(a-r)\frac{\partial}{\partial r}\delta\rho_{2}-\frac{ \omega_{F}^{2}}{4\pi}\delta(a-r)\frac{\partial}{\partial r}\int\frac{d^{2}r_{1 }\delta\rho(\mathbf{r}_{1})}{|\mathbf{r}-\mathbf{r}_{1}|}.\end{array}  \tag{13}\]  
We suppose the solution of the above equation in the form, \(\delta\rho_{2}=\sigma(\Omega,t)\delta(a+0^{+}-r)\) and multiply both sides of this equation by \(r^{2}\) and integrate with respect to \(r\) in arbitrary limits, i.e., \(\int\limits_{l}^{L}r^{2}dr\dots\), such that \(a\in(l,L)\) (this integration removes Dirac deltas), which leads to the equation,  
\[ \begin{array}{l}a^{2}\frac{\partial^{2}\sigma(\Omega,t)}{\partial t^{2}}=- \frac{k_{B}T}{m}\int\limits_{l}^{L}dr\frac{\partial}{\partial r}\left[(n+ \delta\rho_{2})r^{2}\delta(a-r)\right]\\ -\frac{k_{B}T}{m}\sigma(\Omega,t)\int\limits_{l}^{L}r^{2}dr\delta(a-r)\frac{ \partial}{\partial r}\delta(a-r)\\ -\frac{\omega_{F}^{2}}{4\pi}\int\limits_{l}^{L}r^{2}dr\delta(a-r)\frac{ \partial}{\partial r}\int\limits_{0}^{\infty}r_{1}^{2}dr_{1}\int d\Omega\frac{ \delta\rho_{2}(\mathbf{r}_{2})}{|\mathbf{r}-\mathbf{r}_{1}|}\\ -\frac{\omega_{F}^{2}}{4\pi}\int\limits_{l}^{L}r^{2}dr\delta(a-r)\frac{ \partial}{\partial r}\int\limits_{0}^{a}r_{1}^{2}dr_{1}\int d\Omega\frac{ \delta\rho_{1}(\mathbf{r}_{1})}{|\mathbf{r}-\mathbf{r}_{1}|}.\end{array}  \tag{14}\]  
Two first terms in r.h.s. of the above equation vanish, because,  
\[-\frac{k_{B}T}{m}\int\limits_{l}^{L}dr\frac{\partial}{\partial r}\left[(n+ \delta\rho_{2})r^{2}\delta(a-r)\right]=-\frac{k_{B}T}{m}\left[(n+\delta\rho_{2 })r^{2}\delta(r-a)\right]|_{l}^{L}=0 \tag{15}\]and  
\[ \begin{array}{l}\int\limits_{l}^{L}r^{2}dr\delta(a-r)\frac{\partial}{\partial r} \delta(a-r)=a^{2}\int\limits_{l}^{L}dr\frac{1}{2}\frac{\partial}{\partial r} \delta^{2}(a-r)\\ =\frac{a^{2}}{2}\delta^{2}(a-r)|_{l}^{L}=\frac{a^{2}}{2}lim_{\mu\to 0}\frac{1}{ \pi}\frac{\mu}{\mu^{2}+(a-r)^{2}}\delta(a-r)|_{l}^{L}=0.\end{array}  \tag{10}\]  
Two last terms of r.h.s of Equation (A8) can be transformed using the formula [39], \(\frac{1}{\sqrt{1+z^{2}-2 zcos\gamma}}=\sum\limits_{l=0}^{\infty}P_{l}(cos \gamma)z^{l}\), \(for\;\;z<1\), where \(P_{l}(cos\gamma)=\frac{4\pi}{2l+1}\sum\limits_{m=-l}^{l}Y_{lm}(\Omega)Y_{lm}^ {*}(\Omega)\) are Legendre polynomials. This formula leads to the following one,  
\[