## 3 Model Definition
By direct numerical simulations (TDLDA, i.e., the time dependent local density approximation) [3; 4] it has been proved that the volume-surface excitation mixing gradually disappears in larger clusters (with the number of electrons exceeding ca. 60) [3; 4], which supports accuracy of semiclassical RPA description, within which volume plasmons can be separated from the surface ones (even though the latter can be excited by the former ones, due to the last term in Equation (19)). The volume-surface coupling, important for ultra-small radii of metallic clusters, when the shell effects and quantum spill-out are significant, completely disappear, however, with growing sphere dimension to the range of several nanometers in the case metals. We consider here the binary electrolyte system with radius of order of micrometer, when quantum effects specific for ultra-small clusters are negligible. Such an opportunity allows us to formulate an analytical RPA semiclassical description in the form of an oscillator equation, allowing for phenomenological inclusion of damping effects. The energy dissipation effects, though going beyond the RPA model, turned out to be the overwhelming physical property in the case of larger metallic nanospheres [1; 5] (with \(a>10\) nm for Au or Ag) and also for much larger ionic systems, as we will demonstrate it below. Solution of RPA Plasmon Equation: Volume and Surface Plasmon Frequencies in a Spherical Finite Ion System  
Equations (18) and (19) are solved for metallic nanospheres [1] and these solutions can be directly applied to ionic systems. To summarize briefly this analysis, both parts of the plasma fluctuation can represented as follows,  
\[ \begin{array}{l}\delta\dot{p}_{1}(\textbf{r},t)=n\left[f_{1}(r)+F(\textbf{r},t) \right],\ for\ r<a,\\ \delta\dot{p}_{2}(\textbf{r},t)=nf_{2}(r)+\sigma(\Omega,t)\delta(r+\epsilon-a ),\ \epsilon=0+,\ for\ r\geq a,\ (r\to a+),\end{array}  \tag{20}\]  
with initial conditions, \(F(\textbf{r},t)|_{t=0}=0,\ \sigma(\Omega,t)|_{t=0}=0\), (\(\Omega\) is the spherical angle), \((1+f_{1}(r))|_{r=a}=f_{2}(r)|_{r=a}\) (continuity condition), \(F(\textbf{r},t)|_{r=a}=0\), \(\int\rho(\textbf{r},t)d^{3}r=N\) (neutrality condition). Neglecting the static contributions \(f_{1,2}\) describing the spill-out effect, which sharply lowers with the system size growth [1], one can write out the time-dependent parts of the ion concentration fluctuations in the form,  
\[F(\textbf{r},t)=\sum\limits_{l=1}^{\infty}\sum\limits_{m=-l}^{l}\sum\limits_{ i=1}^{\infty}A_{lmn}j_{l}(k_{nl}r)Y_{lm}(\Omega)sin(\omega_{li}t), \tag{21}\]  
and  
\[ \begin{array}{l}\sigma(\Omega,t)=\sum\limits_{l=1}^{\infty}\sum\limits_{m=- l}^{l}\frac{B_{lm}}{a^{2}}Y_{lm}(\Omega)sin(\omega_{0l}t)\\ +\sum\limits_{l=1}^{\infty}\sum\limits_{m=-l}^{l}\sum\limits_{i=1}^{\infty}A_ {lmn}\frac{(l+1)\omega_{p}^{2}}{l\omega_{p}^{2}-(2l+1)\omega_{li}^{2}}Y_{lm}( \Omega)n_{e}\int\limits_{0}^{a}dr_{1}\frac{l^{l+2}_{1}}{a^{l+2}}j_{l}(k_{li}r _{1})sin(\omega_{li}t),\end{array}  \tag{22}\]  
where \(j_{l}(\xi)=\sqrt{\frac{\pi}{2\xi}}l_{l+1/2}(\xi)\) is the spherical Bessel function, \(Y_{lm}(\Omega)\) is the spherical function, \(\omega_{li}=\omega_{p}\sqrt{1+\frac{k_{B}T_{lm}^{2}}{\omega_{p}^{2}a^{2}m}}\) are the frequencies of ion volume self-oscillations (volume plasmon frequencies), \(x_{li}\) are the nodes of the Bessel function \(j_{l}(\xi)\) numerated with \(i=1,2,3\dots\) (cf. Figure 2), \(k_{li}=x_{li}/a\), \(\omega_{l0}=\omega_{p}\sqrt{\frac{l}{2l+1}}\) are the frequencies of ion surface self-oscillations (surface plasmon frequencies). The derivation of the self-frequencies for ionic plasmon oscillations is presented with all details in Appendix A. Amplitudes \(A_{lmi}\) and \(B_{lm}\) are arbitrary in the homogeneous problem and can be adjusted to the initial conditions.