\section{Dipole Modes of metallic spherical shell}
The surface plasmon eigenmodes of the $i$th nanoshell  are described by the equation

\begin{align}
\nabla\left[\theta^{(i)}(\mathbf{r})\nabla \varphi_m^{(i)} \right]= s_m^{(i)} \nabla^2 \varphi_m^{(i)}. 
\end{align}
Function  $\theta^{(i)}(\mathbf{r})$ of the nanoshell is equal to 1 inside the metal of the $i$th nanoshell and 0 elsewhere. 
The $m=\pm1$ dipole eigenmode is given by
\begin{align}
\varphi = \sin\theta e^{im\phi}
\begin{cases}
c_1 \frac{r}{R_1},& 0<r<R_1,\\
c_2 \frac{r}{R_1},+c_3\frac{R_1^2}{r^2},& R_1<r<R_2,\\
c_4 \frac{R_1^2}{r^2},& R_2<r.\\
\end{cases}\label{eq:c_dipole}
\end{align}
One finds  the eigenvalue $s^{(i)}_m$ for these eigenmodes as
\begin{align}
s_m^{(i)}  &= \frac{1\pm\frac{1}{3}\sqrt{1+8\eta^3}}{2},
\label{s_m}
\end{align}
where
\begin{equation}
\eta = \frac{R_1}{R_2}~.
\end{equation}
We impose a normalization condition
\begin{align}
\int_{\mathrm{All~Space}} |\nabla \varphi_m^{(i)} |^2 dv =1~.
\end{align}
There are two solutions corresponding to the $\pm$ signs in Eq.\ (\ref{s_m}). We choose one corresponding to the $-$ sign, which describes a lower frequency mode. Substituting it, we find:
\begin{widetext}
\begin{align}
c_1 &=\left[ \frac{3}{16\pi R_1}\frac{(1-s_m^{(i)})(s_m^{(i)}-\frac{1}{3})}{(s_m^{(i)}-\frac{1}{2})}
\right]^{1/2} =\left[ \frac{-1+4\eta^3+\sqrt{1+8\eta^3}}{16 \pi R_1\sqrt{1+8\eta^3}}
\right]^{1/2},\\
c_2 &=\left(\frac{s_m^{(i)}-\frac{2}{3}}{s_m^{(i)}-1}\right)c_1,\\
c_3 &=\frac{1}{3}\left(\frac{-1}{s_m^{(i)}-1}\right)c_1,\\
c_4 &=\frac{1}{3}\left(\frac{-1}{s_m^{(i)}-\frac{1}{3}}\right)c_1.
\end{align}
\end{widetext}
The electric field of the SP mode inside the gain core, i.e., for $r\le R_1$, is 
\begin{align}
\mathbf{E} &= -A\frac{c_{1} e^{ im\phi}}{R_1}(\mathbf e_{r}\sin\theta+\mathbf e_{\theta}\cos\theta+ i\mathbf e_{\phi})\notag\\
&= -A\frac{c_{1} }{R_1}(\mathbf e_{x}+ im\mathbf e_{y}),
\end{align}
where $\mathbf e_r, \mathbf e_\theta, \mathbf e_\phi$ and $\mathbf e_x, \mathbf e_y$ are the corresponding spherical and Cartesian unit vectors, and
\begin{align}
A & =\sqrt{ \frac{4\pi\hbar s_m^{(i)}}{\epsilon_h s_m^{(i)\prime} }}.
\end{align}

\section{SPP of Plasmonic Array}
Consider the problem of two media with permittivities $\epsilon_d$ and $\epsilon_m$.
The eigenmodes of SP's  in the quasistatic limit can be obtained from 
\begin{align}
\nabla\left[\theta(\mathbf{r})\nabla \varphi\right] = s \nabla^2 \varphi, 
\end{align}
with $s={\epsilon_d}/({\epsilon_d-\epsilon_m})$ and $\theta(\mathbf{r})$ is  1 inside metal and zero elsewhere.
The geometric properties of the eigenmodes are characterized by a generalized eigenvalue problem,
\begin{align}
\nabla\left[\theta(\mathbf{r})\nabla \varphi_n\right] = s_n \nabla^2 \varphi_n, \label{eq:field_eq}
\end{align}
where $s_n$ and $\varphi_n$ are the eigenvalues and eigenfunctions, respectively.
The conditions on the boundary $S$ are
\begin{align}
\varphi_n \cdot\mathbf{n}\nabla\varphi_n =0,
\end{align}
where $\mathbf{n}$ is the normal vector to $S$. 
Multiplying Eq.~\eqref{eq:field_eq} by $\varphi_n$ and integrating by parts, we obtain
\begin{align}
s_n &= \langle n|\theta(\mathbf{r})| n\rangle,\\
|n\rangle &= \nabla \varphi_n,\\
\langle m|n\rangle &= \int d^3r \left(\nabla \varphi_m^*\right)\left(\nabla \varphi_n\right),   
\end{align}
The eigenvectors are orthogonal and normalized as,
\begin{align}
\langle m|n\rangle &= \delta_{mn}.   
\end{align}
Consider that there are $N$ metal particles (same metal).  
\begin{align}
&\nabla\left[\theta(\mathbf{r})\nabla \varphi_n\right] = s_n \nabla^2 \varphi_n,\label{eq:n_particle}\\
&\theta\mathbf{r})=\sum_{i=1}^{N}\theta^{(i)}(\mathbf{r}).
\end{align} 
We can solve Eq.~\eqref{eq:n_particle} in terms of the local basis $\varphi_{m}^{(i)}$, defined as 
\begin{align}
&\nabla\left[\theta^{(i)}(\mathbf{r})\nabla \varphi_{m}^{(i)}\right] = s_{m}^{(i)} \nabla^2 \varphi_{m}^{(i)},\label{eq:local_basis}
\end{align} 
where $\varphi_{m}^{(i)}$ are the eigenmodes of the $i$th metal particle. 
We expand the total field $\varphi_n$ in the local basis
\begin{align}
\varphi_n = \sum_{i ,m}C_{nm}^{(i)}\varphi_{m}^{(i)}\label{eq:l_expansion},
\end{align}
where $C_{nm}^{(i)}$ are the coefficients to be determined. 
Substituting Eq.~\eqref{eq:l_expansion} in  Eq.~\eqref{eq:n_particle} and using Eq.~\eqref{eq:local_basis}
we obtained
\begin{widetext}
\begin{align}
\sum_{j}\sum_{i\neq j,m}C_{nm}^{(i)}\nabla\left[\theta_j(\mathbf{r})\nabla \varphi_{m}^{(i)}\right]
+ 
\sum_{i,m}s_{m}^{(i)}C_{nm}^{(i)}\nabla^2\varphi_{m}^{(i)}= s_n \sum_{i,m}C_{nm}^{(i)}\nabla^2\varphi_{m}^{(i)},
\label{eq:l_expression_0}
\end{align}
\end{widetext}
where the first term describes the inter-particle interaction.
To deal with the first term, we expand the field $\nabla\varphi_{m}^{(i)}$  
 over the other local eigenstates $\nabla\varphi_{m'}^{(j)}$ ($i\neq j$). The field $\nabla\varphi_{m}^{(i)}$  
inside the $j$th particle can be expanded over the local basis $\nabla\varphi_{m'}^{(j)}$, since  
the local basis $\nabla\varphi_{m'}^{(j)}$ forms a complete set of field inside the $j$th metal.
\begin{align}
\nabla\varphi_{m}^{(i)}=\sum_{ m^{\prime}}\beta_{ijmm'}\nabla\varphi_{m'}^{(j)}
\label{eq:l_expression_1}.
\end{align}
The interaction can be written as
\begin{align}
\sum_{j,m'}\sum_{i\neq j,m}C_{nm}^{(i)}\beta_{ijmm'}\nabla\left[\theta_j(\mathbf{r})\nabla \varphi_{m'}^{(j)}\right]
\end{align}
or equivalently
\begin{align}
\sum_{j,m'}\sum_{i\neq j,m}C_{nm}^{(i)}\beta_{ijmm'}s_{m'}^{(j)}\nabla^2\varphi_{m'}^{(j)}.
\end{align}
The eigenmode equation becomes
\begin{widetext}
\begin{align}
\sum_{j,m'}\sum_{i\neq j,m}C_{nm}^{(i)}\beta_{ijmm'}s_{m'}^{(j)}\nabla^2\varphi_{m'}^{(j)}
+ 
\sum_{i,m}s_{m}^{(i)}C_{nm}^{(i)}\nabla^2\varphi_{m}^{(i)}= s_n \sum_{i,m}C_{nm}^{(i)}\nabla^2\varphi_{m}^{(i)}.
\end{align}
\end{widetext}
Multiplying $\phi_{m}^{(j)}$ on the both sides,  integrating by parts and utilizing the othogonality of the eigenmodes, we obtain
\begin{align}
\sum_{i\neq k,m'}C_{nm'}^{(i)}\beta_{ikm'm}s_{m}^{(k)}
=  (s_n-s_{m}^{(k)})C_{nm}^{(k)}
\end{align}
More compactly,
\begin{align}
\sum_{i,m'}\beta_{ikm'm}s_{m}^{(k)}C_{nm'}^{(i)}
=  s_nC_{nm}^{(k)}. \label{eq:SPP_array_main}
\end{align}
Equation \eqref{eq:SPP_array_main} defines an eigenproblem, which yields the egienvalue $s_n$ and the coefficients $C_{nm}^{(k)}$.  
The Hermiticity of $\beta_{ijmm'}$, i.e.,
\begin{align}
\beta_{ijmm'}=\beta_{jim'm}
\end{align}
 is protected by the Green's reciprocity.
The interaction coefficeints $\beta_{ijmm'}$ are given by an integral
\begin{widetext}
\begin{align}
\beta_{ijmm'}& =\frac{1}{s_{m'}^{(j)}}\int \theta^{(j)}(\mathbf{r})\left( \nabla \varphi_{m' }^{(j)*}\right)\left(\nabla\varphi_{m}^{(i)}\right) d^3r=\frac{1}{s_{m}^{(i)}}\int \theta^{(i)}(\mathbf{r}) \left(\nabla \varphi_{m' }^{(j)*}\right)\left(\nabla\varphi_{m}^{(i)}\right) d^3r\\
                    & =\int_{\mathrm{All~Space}} \left(\nabla \varphi_{m' }^{(j)*}\right)\left(\nabla\varphi_{m}^{(i)}\right) d^3r.\label{eq:overplap}
\end{align}
\end{widetext}
Because matrix $\beta_{ikm'm}s_{m}^{(k)}$ in Eq.~\eqref{eq:SPP_array_main} is not Hermitian,  the eigenvectors are not orthogonal, and  the normalization condition is not $\sum_{i,m} |C_{nm}^{(i)}|^2=1 $.
In actuality, the normalization condition is given by
\begin{align}
\sum_{i,k,m',m}C_{nm}^{(k)*}\beta_{ikm'm}C_{nm'}^{(i)}
= 1.\label{eq:normalization}
\end{align}
We can convert Eq.~\eqref{eq:SPP_array_main} into a Hermitian eigenproblem by
a substitution, $\tilde{C}_{nm}^{(i)}=\sqrt{s_m^{(i)}}C_{nm}^{(i)}$.
With this substitution, we rewrite Eq.~\eqref{eq:SPP_array_main} as
\begin{align}
\sum_{i,m'}H_{ikm'm}\tilde{C}_{nm'}^{(i)}
=  s_n\tilde{C}_{nm}^{(k)}~, 
\label{eq:SPP_array_her}\\
H_{ikm'm}=\sqrt{s_m^{(i)}s_{m'}^{(j)}}\beta_{ikm'm}~.
\end{align}
Because $H_{ikm'm}$ is Hermitian, the eigenvectors are orthogonal and the normalization is simply
\begin{align}
\sum_{i,m}|\tilde{C}_{nm}^{(i)}|^2
= 1.\label{eq:normalization_her}
\end{align}



\section{Characters of Topological SP Eigenmodes as  Irreducible Representations of the $D_{3h}$ Point Symmetry Group}

The spasing eigenmodes are surface plasmons (SPs) that should be classified corresponding to irreducible representations of the point symmetry group of the lattice unit cell, which is $D_{3h}$  \cite{Landau_Lifshitz_Quantum_Mechanics:1965}. This group has six such representations: four of them are singlet or one-dimensional (1d), which are $A_1^\prime$ (scalar or fully-symmetric), $A_2^\prime$ (a pseudovector in the $z$-direction),  $A_1^{\prime\prime}$ (pseudoscalar), and $A_2^{\prime\prime}$ (a polar vector in the $z$-direction), and two are doublet or 2d-representations: $E^\prime$ (a polar vector in the $xy$-plane) and $E^{\prime\prime}$ (a pseudovector in the $xy$-plane). We seek for a solution where the local electric field is in the plane of the lattice (the $xy$-plane) and rotates in time for a given lattice site and in space from a site to a site. A basis for such a field can be constituted, e.g., by two spherical harmonics, $Y_{1\pm 1}$ corresponding to left- and right-handed rotations. 

An example of the local fields of the two eigenmodes, which are mutually time-reversed and should be described by the same doublet representation, is provided by our solution of the spaser equations -- see the main text of the article. Here we show such two modes: the $\mathbf K$-mode rotating right [Fig.\ \ref{Spaser_Field_Dynamics}] and the $\mathcal T$-reversed $\mathbf K^\prime$-mode rotating left [Fig.\ \ref{Spaser_Field_Dynamics_L}].

To unambiguously identify the corresponding representation, we note that the underlying material system is a sublattice A whose unit cell has the third-oder vertical ($z$) axis ($3C_3$ symmetry elements), three second-order axes (in the $xy$ plane) ($3U_2$ symmetries), also $3S_3$ elements (mirror-rotational symmetries), and $3\sigma_v^\prime$ elements (reflections in the vertical planes). 
Spherical harmonics $Y_{1\pm1}$ are spherical representations of a polar vector in the $xy$ plane: 
\begin{equation}
Y_{1\pm1}=\mp\frac{1}{2}\sqrt{\frac{3}{2\pi}}(x\pm iy)~.
\label{Y1pm1}
\end{equation}

\begin{figure}
\includegraphics[width=0.5\columnwidth]{Spaser_Field_Dynamics.png}
\caption{Dynamics of $\mathbf K$-valley SP eigenmode. (a)-(d) The distributions of the local field modulus are plotted for the unit cell in the real space. The phases of the spaser oscillation are indicated for each panel. The color-coded scale of the local fields is given at the left=hand side for one SP quantum per the unit supercell. The local fields rotate in time  in the direction of the black arrow, i.e., clockwise. (e) Edge fields for $\mathbf K$-valley mode. The instantaneous dipole vectors are indicated by the green arrows; the direction of the field rotation is depicted by the curved black arrows;  crystal momentum $\mathbf k$ of the mode is shown by the straight black arrow. Each such a rotating dipole generates a current that is normal to this dipole.
\label{Spaser_Field_Dynamics}
}
\end{figure}
\begin{figure}
\includegraphics[width=0.5\columnwidth]{Spaser_Field_Dynamics_L.png}
\caption{Dynamics of the $\mathbf K^\prime$-valley SP eigenmode. Same as Fig.\ \ref{Spaser_Field_Dynamics} but the dynamics is $\mathcal T$-reversed: fields rotate left and propagate left.
\label{Spaser_Field_Dynamics_L}
}
\end{figure}
\begin{table}[]
\begin{tabular}{@{}ccccccc@{}}
\hline
                              & $E$                   & $\sigma_h$ & $3C_3$  & $3S_3$  & $3U_2$  &  $3\sigma_v^{\prime}$  \\ \hline
$E^{\prime}$ &    2                   &  2              & -1         & -1           &  0          &    0   \\ \hline
\end{tabular}
\caption{Characters of the $E^{\prime}$ representation of the $D_{3h}$ point symmetry group, adapted from Ref. \onlinecite{Landau_Lifshitz_Quantum_Mechanics:1965}.
The elements of $D_{3h}$ are the identity $E$, the horizontal-plane reflection $\sigma_h$, the three rotations $C_3$, the three rotary-reflections $S_3$, 
the three rotations about the horizontal axes $U_2$, and the three vertical-plane reflections $\sigma_v^{\prime}$.
}\label{tb:ch}
\end{table}

The character of the trivial $E$ operation is always the representation dimensionality, i.e., $\chi(E)=2$. To find $\chi(\sigma_h)$ for our eigenvectors we note that the $\sigma_h$ operation is reflection in the $xy$ mirror plane (i.e., a transformation $x\to x$, $y\to y$, and $z\to -z$), which leaves both of the basis functions of Eq.\ (\ref{Y1pm1}) unchanged. Correspondingly, $\chi(\sigma_h)=2$. Operation $C_3$ is a rotation by a $2\pi/3$ angle whose character, obviously, is $\chi(C_3)=2\cos\left(2\pi/3\right)=-1$. An operation $S_3$ is a product of $\sigma_h$, which leaves the basis unchanged, and $C_3$; thus its character is $\chi(S_3)=\chi(C_3)=-1$. A $U_2$ rotation, say about the $y$ axis of symmetry, causes $y\to y, x\to -x$; this results in $Y_{11}\leftrightarrow Y_{1-1}$. Thus the basis functions are exchanged, and the corresponding transformation matrix is purely off-diagonal; correspondingly, $\chi(U_2)=0$. An operation $\sigma_v^\prime$ acts for our purposes exactly like $U_2$; thus $\chi(\sigma_v^\prime)=0$. These characters  that we have obtained are exactly those of the $E^\prime$ representation of the $D_{3h}$ group summarized in Table \ref{tb:ch}. This unambiguously identifies that our spasing eigenmodes transform accordingly to the $E^\prime$ doublet representation of the $D_{3h}$ point symmetry group. Note that that the two components of this doublet describe field distribution rotating in time and space in the opposite directions. They transform one into another by either of the following transformations: $\mathcal T$, $\sigma_v^\prime$, or $U_2$. 


