\section{Introduction}
Strong coupling between two harmonic oscillators -- either of classical or quantum nature -- is one of the most basic physical models that can be employed to understand the behavior of various mechanical and electromagnetic systems.
Polaritons -- stationary eigenstates of a coupled system in the strong coupling regime -- are hybridized states whose wave function is characterized by the photonic and the matter component simultaneously \cite{Mills1974,torma2014strong}.
In the optical domain, polaritonic states are often realized by means of coupling an optical cavity mode with excitonic or vibrational transitions in resonant media \cite{khitrova2006vacuum, torma2014strong, Baranov2018}.
Thanks to their hybrid composition, optical polaritons manifest unique properties that are typical of simultaneous excitations of light and matter \cite{sanvitto2016road}. This offers new ways for modifying the microscopic properties of matter \cite{Galego2015,ebbesen2016hybrid} and even controlling the rates of chemical reactions \cite{thomas2016ground,Herrera2016,Munkhbat2018,Thomas2019,peters2019effect,stranius2018selective}. The extent to which these microscopic properties can be modified is often determined by the interaction strength between the two components of the coupled system \cite{galego2016suppressing,martinez2018can, feist2018polaritonic,fregoni2022theoretical,schafer2019modification}. This calls for finding the ways to improve the coupling strength in polaritonic systems.



In the saturation limit -- when the optical cavity is saturated with resonant transitions -- the key characteristic that determines the resulting interaction constant is the reduced oscillator strength, which is proportional to the transition dipole moment, and the volume density of oscillators in the medium \cite{Platts2009,canales2021abundance}. There is a solid theoretical evidence to believe that it is not possible to boost the coupling strength by confining the cavity mode volume -- shrinking the cavity mode volume increases the coupling strength with few to individual emitters \cite{chikkaraddy2016single}, up to the point where the mode will excite less and less transitions as the volume becomes smaller than the emitters themselves \cite{rossi2019strong,Kuisma2022}.


Recently, it was proposed to boost/increase the coupling strength in polaritonic systems by utilizing so called meta-atoms -- resonant metallic nanoparticles hosting localized plasmonic resonances \cite{ameling2010cavity, bisht2018collective}. Despite not having a discrete anharmonic energy ladder like quantum emitters, such meta-atoms participate in the coupling process in a similar way, resulting in the emergence of hybrid polaritonic states with equally spaced (harmonic) energy ladders \cite{konrad2015strong, hertzog2021enhancing}. This approach has enabled the exotic regimes of ultra-strong (USC) \cite{Baranov2020,rajabali2022ultrastrongly} and even deep strong coupling \cite{mueller2020deep} at ambient conditions, which was previously unavailable with more traditional quantum emitter platforms. In these regimes, not only the excited states, but also the ground state of the system experiences a modification upon coupling \cite{ciuti2005quantum,de2014light,forn2019ultrastrong}.


The experimental progress in realizing polaritonic states with artificial meta-atoms begs a natural question: what is the upper bound on the coupling strength in polaritonic systems involving meta-atoms?
Some analytical models have been proposed that describe polaritonic states of such system using rigorous coupled-dipoles method \cite{lamowski2018plasmon}. However, that particular model was developed for subwavelength spherical meta-atoms, and as a result cannot be applied to study the limits of plasmon-photon coupling in the case of large spherical or non-spherical particles.

\begin{figure}
\centering\includegraphics[width=1\columnwidth]{Fig_1.pdf}
\caption{Geometry of the system under study: an ensemble of spheroidal metallic nanoparticles (meta-atoms) described by Drude permittivity $\varepsilon_m$ distributed in the background medium (air) with a volume filling factor $f$. The ensemble of meta-atoms interacts with the photonic mode propagating through the background medium (air).}
\label{fig1}
\end{figure}

In this paper, using a Hamiltonian approach we theoretically analyze the formation of polaritonic states and examine the upper limits of the collective plasmon-photon coupling strength in a number of dense assemblies of plasmonic meta-atoms. Starting with the case of analytically solvable spherical meta-atoms, we identify the universal upper bounds on the normalized collective coupling strength $g/\omega_0$ between ensembles of plasmonic meta-atoms and free-space photons. Then, with the aid of numerical simulations, we examine the case of spheroidal meta-atoms. Our results suggest that strongly elongated meta-atoms could be optimal geometry for attaining the highest value of the collective coupling strength with the optical field.




\section{System under study}


The system we analyze in this study is represented by an ensemble of (generally) spheroidal metallic meta-atoms with semi-axes $a$, $b$, and $c$, and permittivity $\varepsilon_m(\omega)$ distributed in a host medium (air), Fig. \ref{fig1}.
The meta-atoms are assumed to be distributed in space  periodically, with volume filling factor $f$ defined as a ratio of the volume occupied by the metallic fraction to the total volume of the system, $f=V_{metal}/V$.


The permittivity of metallic meta-atoms is described by the Drude model:
\begin{equation}
    \varepsilon_m (\omega) = 1 - \frac{\omega_p^2}{\omega (\omega + i \gamma)},
    \label{Eq_2}
\end{equation}
where $\omega_p$ is the plasma frequency of the metal. For the sake of simplicity, we will ignore dissipation in this study and assume $\gamma = 0$. 
Each meta-atom supports a localized plasmon resonance associated with oscillations of the electron density of the particle.
The resulting ensemble of plasmonic meta-atoms interacts with a photon of energy $\hbar \omega_{\mathbf{k} }$ propagating across the surrounding medium, which is assumed to be air ($\varepsilon = 1$).

In order to describe the polaritonic spectrum of the system, we will develop in the following a Hamiltonian model involving the matter polarization. For simplicity, we describe the interaction of the localized oscillators' dipolar transitions with the transverse photonic field, which can be generalized to multipolar transitions.

\section{Hamiltonian model of the system}

\begin{figure*}
    \centering
    \includegraphics[width=1\textwidth]{Fig_2.pdf}
    \caption{The general Power-Zienau-Woolley Hamiltonian (Eq. \eqref{Hamfull}) is used to describe a single meta-atom with polarization $\mathbf{P}(\mathbf{r},t)=\alpha{\bf X}(\mathbf{r},t)$. Assuming negligible Joule losses, the dipole moment $\boldsymbol{\mu}_{0,1}$ is obtained with Fermi's golden rule via the linewidth $\gamma_{pl}$. The same general Hamiltonian is then used in the continuous limit to model the meta-atom array as a density of dipoles per unit volume $V_\text{cell}$, with inputs $\omega_{pl}$ and $\boldsymbol{\mu}_{0,1}$ from the single meta-atom description.}
    \label{H_meta_atoms}
\end{figure*}

  \subsection{General Power-Zienau-Woolley Hamiltonian}
  
To describe the interaction between the cubic lattice of meta-atoms and the transverse photonic field, we use a Hamiltonian formulation in the Power-Zienau-Woolley representation of the Coulomb gauge. The authors of ref. \cite{lamowski2018plasmon} previously developed a  model of the plasmon-photon interaction in the Coulomb gauge that accurately takes into account the near-field Coulomb interaction between meta-atoms. The resulting model, in particular, captures the anisotropy of the polaritonic modes induced by the internal anisotropy of the cubic lattice of the meta-atoms. However, the limits to the plasmon-photon coupling have not been analyzed.

Here, we are aiming at a simplified model that would allow us to easily characterize the strength of plasmon-photon coupling in densely packed cubic arrays. To that end, we employ the Power-Zienau-Woolley (or multipolar) representation, wherein the dipole-dipole interactions between the meta-atoms are accounted for by the quadratic self-polarization term. Considering an ensemble of meta-atoms, each of volume $V_j$ and single resonance frequency $\omega_j$, we show in the Supplementary Information that the Hamiltonian:
\begin{multline}
\label{Hamfull}
H = \underbrace{\sum_j\int_{V_j}\mathrm{d}^3r\left(\frac{\boldsymbol{\Pi}_j^2}{2\rho} + \frac{1}{2}\rho\omega_{j}^2{\bf X}_j^2\right)}_\text{meta-atoms} \\+ \underbrace{\int_{\mathbb{R}^3}\mathrm{d}^3r\left(\frac{\boldsymbol{\Pi}_A^2}{2\epsilon_0} + \frac{(\nabla \times \mathbf{A})^2}{2\mu_0}\right)}_\text{radiation}\\ + \underbrace{\int_{V_m}\mathrm{d}^3r\frac{\mathbf{P}^2}{2\epsilon_0}}_\text{$\mathbf{P}^2$-term} + \underbrace{\int_{V_m}\mathrm{d}^3r\frac{\mathbf{P}\cdot\boldsymbol{\Pi}_A}{\epsilon_0}}_\text{interaction}
\end{multline}
generates the source-free Maxwell's equations $\nabla\cdot\mathbf{D}=0,\nabla\cdot\mathbf{B}=0,\nabla \times \mathbf{E}=-\frac{\partial}{\partial t}\mathbf{B},\nabla \times \mathbf{H}=\frac{\partial}{\partial t}\mathbf{D}$, as well as the equation of motion of the macroscopic polarization $\mathbf{P}=\sum_j\mathbf{P}_j$:
\begin{align}
    \left(\frac{\partial^2}{\partial t^2} + \omega_{j}^2\right)\mathbf{P}_j(\mathbf{r},t) = \beta_j(\mathbf{r})\mathbf{E}(\mathbf{r},t),
\end{align}
where $\beta_j(\mathbf{r}) = \alpha_j(\mathbf{r})/\rho$, $\rho$ being the carrier volumic mass in the medium of the resonators, and $\alpha_j(\mathbf{r})$ being the displaced carrier density in resonator $j$. The displacement of the carriers is contained in the vector field ${\bf X}_j(\mathbf{r},t)$, whose canonical momentum is $\boldsymbol{\Pi}_j=\rho\frac{\partial}{\partial t}{\bf X}_j$. The radiation field with vector potential $\mathbf{A}$ and canonical momentum $\boldsymbol{\Pi}_A=-\mathbf{D}$ is purely transverse in the Coulomb gauge. We emphasize that, in the picture given by Hamiltonian \eqref{Hamfull}, Ohmic dissipation is disregarded. A more rigorous treatment with additional degrees of freedom in the Hamiltonian would provide a description of the Ohmic losses, but is beyond the scope of this work. The general form of Hamiltonian \eqref{Hamfull} is used to treat both the cases of: 1) a periodic array of densely-packed meta-atoms, and 2) the description of a single meta-atom coupling to the radiation field, enabling us to input parameters from the single meta-atom in the collective description.

  \subsection{Densely packed meta-atom array Hamiltonian}

Assuming a densely packed array, we perform the continuous limit, i.e. we describe the array as a density of dipoles with a polarization expanded in the plane wave basis \cite{todorov2014dipolar}:
\begin{align}
\mathbf{P}(\mathbf{r}) = \sum_\mathbf{Q}\frac{\boldsymbol{\mu}_N}{V_\text{cell}}\mathrm{e}^{-i\mathbf{Q}\cdot\mathbf{r}}\left(b_\mathbf{Q} + b_{-\mathbf{Q}}^\+\right),
\end{align}
with $\boldsymbol{\mu}_N = \sqrt{N}\boldsymbol{\mu}_{0,1}$ being the effective dipole moment in the lattice volume $V_\text{cell}$, containing $N= fV_\text{cell}/(\frac{4}{3}\pi r^3)$ meta-atoms. In the above formula, $b_\mathbf{Q},b_\mathbf{Q}^\+$ are bosonic annihilation and creation operators for a matter excitation propagating with wavevector $\mathbf{Q}$ in the array. The matter Hamiltonian $H_\text{mat}$, containing the bare meta-atom and $\mathbf{P}^2$ parts of \eqref{Hamfull} is then, for a given $\mathbf{Q}$:
\begin{align}
    H_\text{mat} =\hbar\omega_{pl}b_\mathbf{Q}^\+b_\mathbf{Q} + \frac{\boldsymbol{\mu}_N^2}{2\epsilon_0V_\text{cell}}\left(b_\mathbf{Q} + b_{-\mathbf{Q}}^\dag\right)\left(b_{-\mathbf{Q}} + b_\mathbf{Q}^\dag\right).
\end{align}
The diagonalization of the matter Hamiltonian $H_\text{mat} \rightarrow \hbar\Omega_{pl}B_\mathbf{Q}^\dag B_\mathbf{Q}$, yields the collective eigenfrequency of the meta-atom array: $\Omega_{pl} = \sqrt{\omega_{pl}^2+2\boldsymbol{\mu}_N^2\omega_{pl}/(\hbar\epsilon_0V_\text{cell})}$.

We next focus on the light-matter interaction term, assuming the array is homogeneous in the dipole orientation, along the $z$-axis. In the original basis, involving the operators $b_\mathbf{Q},b_\mathbf{Q}^\+$, the light-matter interaction term corresponding to the last term of Eq. \eqref{Hamfull} takes the form $H_\text{light-mat} = -i\hbar g_\mathbf{k}  (\widehat{a} _{\mathbf{k} } - \widehat{a} _{-\mathbf{k} }^{\dag})(\widehat{b} _{-\mathbf{k} } + \widehat{b} _{\mathbf{k} }^{\dag})$, where $a_{\mathbf{k} }$ is the annihilation operator for a transverse magnetic photon with wave vector $\mathbf{k} $ and frequency $\omega_{k} = c|\mathbf{k} |$, and the array excitation wavevector $\mathbf{Q}$ must now match the photon wavevector $\mathbf{k}$. The \emph{collective} coupling constant is $\hbar g_{\mathbf{k} } = \mathcal{E}_\text{vac}(\mathbf{k} )\boldsymbol{\mu}_N\cdot\boldsymbol{\epsilon}_\mathbf{k} $, where $\mathcal{E}_\text{vac}(\mathbf{k} ) = \sqrt{\hbar \omega_{k}/ (2\varepsilon_0 V_\text{cell})}$ is the vacuum electric field of the photonic mode confined to the quantization box of volume $V_\text{cell}$ and $\boldsymbol{\epsilon}_\mathbf{k} $ is the unit transverse magnetic polarization vector.
One can see that the volume $V_\text{cell}$ is cancelled out from the resulting expression for the coupling strength $g_\mathbf{k} $, which now depends on the meta-atom's transition dipole moment and the meta-atoms density: 
\begin{equation}
    g_{\mathbf{k} } = \sqrt{\frac{f}{V_0} \frac{\hbar \omega_{k}}{2\varepsilon_0}}\boldsymbol{\mu}_{0,1}\cdot\boldsymbol{\epsilon}_\mathbf{k} 
    \label{Eq_9}
\end{equation}
where $V_0$ is the physical volume of a single meta-atom. In the new basis involving the diagonalized matter part, the Hamiltonian is, for a given wavevector $\mathbf{k}$:
\begin{multline}
\label{H_TM}
    H = \hbar\omega_ka_\mathbf{k}^\+a_\mathbf{k} + \hbar\Omega_{pl}B_\mathbf{k}^\+B_\mathbf{k} \\- i\widetilde{g}_\mathbf{k}\left(a_\mathbf{k} - a_{-\mathbf{k}}^\+\right)\left(B_{-\mathbf{k}} + B_{\mathbf{k}}\right),
\end{multline}
where $\widetilde{g}_\mathbf{k}  = g_\mathbf{k} \sqrt{\omega_{pl}/\Omega_{pl}}$ is the rescaled coupling strength resulting from the first diagonalization step. The full diagonalization of the light-matter Hamiltonian is detailed in the Supplementary Information. It is based on determining the eigenvalues $\Omega_\pm(\mathbf{k})$ of the Hopfield matrix associated to the Hamiltonian \eqref{H_TM}, which are then:
\begin{align}
    \Omega_\pm(\mathbf{k} ) = \frac{1}{\sqrt{2}}\sqrt{\omega_k^2+\Omega_{pl}^2\pm\sqrt{\left(\omega_k^2-\Omega_{pl}^2\right)^2+16\widetilde{g}_\mathbf{k} ^2\omega_k\Omega_{pl}}}.
    \label{Eq_polaritons}
\end{align}


\subsection{Transition dipole moments of individual meta-atoms}

One important ingredient in the developed model is the transition dipole moment of meta-atoms. To that end, we first consider a single meta-atom in free space and employ the classical theory of light scattering in order to evaluate its transition dipole moment matrix element.

The localized plasmon quasi-normal mode of a single plasmonic meta-atom can essentially be described as a harmonic oscillator with an equidistant energy ladder (see Fig. \ref{H_meta_atoms}).
In the spirit of cavity QED dealing with subwavelength two-level quantum emitters, such as electronic and vibrational transitions of atoms and molecules, we will characterize the transitions between each pair of eigenstates $n \to n+1$ by the transition dipole moment (TDM) matrix element $\boldsymbol{\mu}_{mn}$. 
We restrict our analysis to the low-energy domain of the oscillator represented by the $0 \to 1$ transition of the meta-atoms.
Particularly, the radiative transition between the ground $\rvert 0 \rangle$ and the first excited $\rvert 1 \rangle$ Fock states of the plasmonic meta-atom is quantified by:
\begin{equation}
    \boldsymbol{\mu}_{0,1} = \langle 0 \rvert q \mathbf{r}  \rvert 1 \rangle,
\end{equation}
where $\mathbf{r} $ is the position operator corresponding to the center-of-mass motion of the electron cloud for the dipolar plasmon, with total charge $q$.


\begin{figure}[b!]
\centering\includegraphics[width=1\columnwidth]{Fig_3.pdf}
\caption{(a) Magnitude of the transition dipole moment matrix element $|\boldsymbol{\mu}_{0,1}|$ of spherical plasmonic meta-atoms as a function of the meta-atom radius $r$ described by the Drude permittivity, Eq. \ref{Eq_2}, for $\omega_p=5$ eV calculated as poles of the characteristic equation, Eq. \ref{Eq_3}. (b) Normalized transition dipole moments $|\boldsymbol{\mu}_{0,1}|/r$ as a function of the dimensionless radius $r/\lambda_p$ for three values of the plasma frequency.}
\label{fig2}
\end{figure}



Although the system, generally, consists of spheroidal meta-atoms, we begin our analysis with the special case of spherical meta-atoms, whose properties can be described by closed-form analytical expressions.
To quantify these dipole transitions, we first find complex eigenfrequencies of the TM$_1$ quasinormal modes of the metallic sphere.
These eigenfrequencies can be found numerically as roots of the characteristic equation \cite{Bohren2004}:
\begin{equation}
    {n \psi_l(n x)  \xi_l '(x) - \xi_l (x)\psi_l '(n x)}=0,
    \label{Eq_3}
\end{equation}
where $x=k_0 r$, $\psi_l(x)=x j_l(x)$ and $\xi_l(x)=x h_l^{(1)}(x)$ are Ricatti-Bessel functions, and $j_l(x)$ and $h_l^{(1)}(x)$ are spherical Bessel and Hankel functions of the first kind, respectively.


\begin{figure}[b!]
\centering\includegraphics[width=.5\columnwidth]{Fig_4.pdf}
\caption{Geometry of the periodic system incorporating spherical meta-atoms: a periodic ensemble of metallic spherical meta-atoms of radius $r$ forming an FCC lattice with the lattice constant $a$.}
\label{fig4lattice}
\end{figure}


Once the complex-valued eigenfrequencies $\widetilde{\omega}=\omega_0 - i \gamma/2$ of the electric dipole quasinormal modes of the meta-atoms are determined, one can calculate their transition dipole moments $\mu=|\boldsymbol{\mu}_{0,1}|$ by applying the spontaneous decay rate formula \cite{novotny2012principles}:
\begin{equation}
    \gamma = \frac{\omega_0^3}{3 \pi \hbar \varepsilon_0 c^3} |\boldsymbol{\mu}_{01}|^2.
    \label{Eq_4}
\end{equation}
Although this expression is traditionally used to describe radiative decay rate of two-level systems, it can be equally applied to describe the transition rates between the equidistant levels of a harmonic multi-level emitter. Using Wigner-Weisskopf theory, we show in the Supplementary Information that, restricted to a single meta-atom, the Hamiltonian \eqref{Hamfull} yields Eq. \eqref{Eq_4}.


Fig. \ref{fig2}(a) shows the resulting transition dipole moments of a spherical metallic meta-atom as a function of the meta-atom radius evaluated for the value of plasma frequency $\omega_P = 5$ eV. 
The plot reveals a nearly linear dependence of the transition dipole moment on the meta-atom radius. The result for a series of other plasma frequencies (10 and 15 eV) shows an analogous behavior (see Fig. S2). This linear dependence is confirmed by the fit of the dipole moment in the double-logarithmic scale, Fig. S2(c). We have checked that the linear dependence remains valid at large meta-atom radius, up to 5 microns, see Fig. S3. 
For completeness, we also show in Fig. S4 the corresponding resonant frequencies $\Re[\widetilde{\omega}] = \omega_0$ of the electric dipole quasinormal modes of meta-atoms.

\begin{figure*}[hbt!]
\centering\includegraphics[width=.8\textwidth]{Fig_5.pdf}
\caption{(a) Collective light-matter coupling strength $g_{\mathbf{k} }$ in the array of densely packed spherical plasmonic meta-atoms as a function of the meta-atom radius $r$. (b) Normalized collective coupling strength $\eta = g_{\mathbf{k} }/\omega_0$ for spherical meta-atoms as a function of the dimensionless radius $r/\lambda_p$. The values for a series of plasma frequencies are shown. (c) An exemplary spectrum of polaritonic eigenenergies, Eq. \ref{Eq_polaritons}, calculated for an array of spherical meta-atoms for $\omega_p = 10$ eV and the meta-atom radius $r=500$ nm. Shaded area denotes the polariton gap with no allowed real-valued energies. (d) Normalized value of the polariton gap $\Delta/\omega_0$  as a function of the dimensionless radius $r/\lambda_p$ for three values of the plasma frequency.}
\label{fig3}
\end{figure*}

One can notice that for any fixed value of filling factor $f$ there are only two independent dimensional parameters with units of length completely determining the problem with an individual meta-atom: the meta-atom's radius $r$, and the plasma frequency $\omega_P$, which can be translated to the corresponding plasma wavelength $\lambda_P = 2\pi c/\omega_P$. This suggests that the normalized radius $r/\lambda_P$ may play the role of a dimensionless parameter, which could determine all dimensionless characteristics of the meta-atom. To verify this hypothesis, we plot in Fig. \ref{fig2}(b) the normalized transition dipole moment $\mu/r$ as a function of the dimensionless radius $r/\lambda_P$, and merge the data series for all three plasma frequencies we have studied. The result clearly shows that the normalized transition dipole moments follow the same  dependence with variation of $r/\lambda_P$, and approach a universal constant for large $r/\lambda_P$ irrespective of the plasma frequency of the underlying material.



Overall, the data presented above strongly indicates that the transition dipole moment of large plasmonic meta-atoms grows linearly with radius. To substantiate this asymptotic behavior, we propose a simple analytical estimation of the transition dipole moment.
In the limit of large radii of the metallic sphere its fundamental TM1 resonance gradually shifts towards long wavelengths, $\lambda_0 \to \infty$. 
In this limit Drude metal described by Eq. \ref{Eq_2} turns into a perfect electric conductor with $\Re[\varepsilon] < 0,\ |\varepsilon| \gg 1$, and the metallic nanoresonator turns into a dipole antenna \cite{barnard2008spectral}. 
Correspondingly, all dimensional characteristics of the antenna's resonance - such as the resonant wavelength and inverse linewidth - begin to scale linearly with the radius \cite{novotny2007effective, milligan2005modern}: $\lambda_0 = A r$, $2\pi c/\gamma = B r$, where $A$ and $B$ are dimensionless parameters.





Since they are dimensionless and at this point there is only one quantity with the units of length left -- plasma wavelength $\lambda_P$ -- the constants $A$ and $B$ cannot depend on the plasma frequency of the Drude metal.
In turn, this scaling leads to a constant quality factor of the resonance for large radii:
\begin{equation}
    Q = \frac{\omega_0}{\gamma} \to  \frac{B}{A} = \textrm{const},
\end{equation}
which is independent of the plasma frequency. Figure S5, showing the radius dependence of the $Q$-factors of spherical meta-atoms, confirms this statement.


Plugging these asymptotic dependencies into Eq. \ref{Eq_4} and resolving it with respect to $\mu$, we obtain the linear scaling of transition dipole moment of the metallic sphere with radius:
\begin{equation}
    \mu = \sqrt{\frac{A^3}{B}\frac{3 \hbar \varepsilon_0 c}{4 \pi}}  r,
\end{equation}
which agrees with the behavior shown in Fig. \ref{fig2}.
Alternatively, we can express the transition dipole moment from Eq. \ref{Eq_4} in terms of the meta-atom $Q$-factor. The result takes the form:
\begin{equation}
    \mu = \sqrt{AB\frac{3 \hbar \varepsilon_0 c}{4 \pi}} \frac{r}{Q}.
\end{equation}
This expression will be useful in the following analysis of the collective coupling strength in the arrays of meta-atoms.
Additionally, Fig. S6 shows the diameter-to-resonant wavelength ratio $2r/\lambda_0$ for spherical meta-atoms, confirming linear scaling of the resonant wavelength with radius.









\section{Results}

\subsection{Polaritonic spectra with spherical meta-atoms}

Having calculated transition dipole moments of individual meta-atoms, we now turn to the analysis of the collective polaritonic states in densely packed meta-atom ensembles. 
For the following we assume that spherical meta-atoms occupy the sites of a three-dimensional face centered cubic (FCC) lattice with a lattice constant $a$, Fig \ref{fig4lattice}.  For the FCC lattice of spheres of radius $r$ and center-to-center distance $a$, the filling factor is given by
\begin{equation}
    f =  4 \dfrac{4 \pi r^3}{3a^3}.
    \label{Eq_1_f}
\end{equation}
The closest configuration of spheres in such lattice can be reached for $a=2r\sqrt{2}$, when Eq. \ref{Eq_1_f} recovers the famous filling factor of an array of closely packed spheres, $f_0 = \pi/(3\sqrt{2}) \approx 0.74$.
The resulting system closely resembles the structure realized in ref. \cite{mueller2020deep}.


Fig. \ref{fig3}(a) shows the resulting unscaled collective plasmon-photon coupling strength $g_{\mathbf{k} }$ as a function of the meta-atom radius evaluated for $\omega_p=10$ eV at the zero-detuning condition ($\omega_{ph} = \omega_{pl}$). The resulting coupling strength monotonically decreases with the radius of single meta-atom. A more interesting behavior is found when we plot the \emph{normalized} coupling strength $\eta=g_{\mathbf{k} }/\omega_{pl}$ versus the normalized radius $r/ \lambda_P$ for all three plasma frequencies, Fig. \ref{fig3}(b). 


\begin{figure}[hbt!]
\centering\includegraphics[width=.8\columnwidth]{Fig_6.pdf}
\caption{Mode splitting $\Omega$ in an array of $r=500 nm$ and $\omega_p=10$ eV spherical meta-atoms as a function of the center-to-center distance $a/r$. Even for diluted structures with $a/r > 10$ the mode splitting remains well above the characteristic non-radiative plasmon decay rate of 50 meV.}
\label{fig3b}
\end{figure}


As in Fig. \ref{fig2}(b), normalized coupling strengths obtained for different plasma frequencies follow the same common dependence. 
This plot clearly shows that even relatively small meta-atoms easily reach the regime of ultrastrong light-matter coupling with $\eta > 0.5$.  This is the regime of interaction in which the standard quantum optical approximations, such as the rotating wave approximation, fail. Thus so-called fast-rotating terms, as well as the quadratic $P^2$ term (or, alternatively, the $A^2$ term of the Hamiltonian before the PZW transformation) must be taken into account in order to correctly describe the system's behavior \cite{ciuti2005quantum,Todorov2012,todorov2015dipolar,schafer2020relevance}. 
More remarkably, Fig. \ref{fig3}(b) reveals even more exotic domain of deep strong coupling, commonly defined as the regime of interaction with $\eta>1$ \cite{DeLiberato2014, casanova2010deep, langford2017experimentally}. 
From the obtained data we can conclude that the transition to this regime occurs close to $r = \lambda_P$; whether this is the exact threshold or not, should be the subject of a more accurate analytical treatment.


\begin{figure}[b!]
\centering\includegraphics[width=.9\columnwidth]{Fig_7.pdf}
\caption{The attainable maximum $\Bar{\eta}$ of normalized collective coupling strength between the fundamental electric dipole transition of spheroid meta-atoms and free-space photonic field for closely packed arrays of metallic spheroids  as a function of the spheroid AR. The data is obtained for a fixed plasma frequency $\omega_p = 10$ eV. Inset: geometry of the modified meta-atom, represented by a prolate spheroid metallic meta-atoms with the long axis $a$ and shorter axis $b$.}
\label{fig4}
\end{figure}


Another remarkable feature of the data in Fig. \ref{fig3}(b) is that the normalized coupling strength $\eta$ asymptotically approaches a constant (the same for all plasma frequencies) in the limit of large radius:
\begin{equation}
    \eta \le 1.2.
\end{equation}
Since the data points obtained for different plasma frequencies follow the same dependence, this suggests that this upper bound is universal for all plasma frequencies and depends only on the filling factor and the meta-atom shape, which we are going to address below.


This asymptotic behavior can be understood on account of the linear scaling of transition dipole moments in the limit of large radii that is reported in Fig. \ref{fig2}.
Indeed, taking into account the asymptotic behaviors of the resonant meta-atom energy $\omega_0$, transition dipole moment $\mu$, and combining it with the physical meta-atom volume $V_0$, we obtain that the normalized coupling strength approaches a constant:
\begin{equation}
    \frac{g_{\mathbf{k} }}{\omega_0} \to \frac{3 A \hbar}{8 \pi} \sqrt{\frac{B}{\pi}} \frac{\sqrt{f}}{Q} = \textrm{const}.
     \label{Eq_17}
\end{equation}
This simple argument does not, however, allow us to estimate the exact value of the upper limit for the normalized collective coupling strength.



Using the obtained coupling strengths we present a typical spectrum of polaritonic eigenenergies $\Omega_{\pm}$ of an fcc array of meta-atoms with $\omega_P=10$ eV and $r=500$ nm, Fig. \ref{fig3}(c). The dispersion features a familiar anti-crossing picture with a mode splitting of $\Omega_R = 2g \approx 0.8$ eV.
Fig. \ref{fig3}(c) also exhibits a polaritonic gap -- a region of energies with no polaritonic states within it \cite{hopfield1958theory}.
The polariton gap  can be interpreted as the Reststrahlen band of the material: the domain of energies wherein the real part of the permittivity becomes negative, thus forbidding propagation of plane waves \cite{canales2021abundance,Todorov2012}. The lower edge of this gap is exactly the uncoupled meta-atom's energy $\omega_{pl}$. The upper edge of the polariton gap is obtained by calculating the upper polariton energy in the limit $k = 0$ and is $\omega = \sqrt{\omega_0^2 + 4g^2}$. The width of the polariton gap therefore is
\begin{equation}
    \Delta = \sqrt{\omega_0^2 + 4g^2} - \omega_0.
\end{equation}
Fig. \ref{fig3}(d) presents the normalized width of the polariton gap $\Delta_{pol}/\omega_{pl}$ as a function of the normalized radius $r/ \lambda_P$ for the three studied plasma frequencies. Like in other instances, all data series follow a common dependence, once again highlighting the key role of the dimensionless radius of the meta-atom in this coupling problem.

Although our Hamiltonian model does account for the inter-particle Coulomb interactions, it may become less suitable for closely packed lattices with touching particles, when the higher-order multipole interactions beyond the dipole-dipole one may become dominating. 
For this reason we now apply it for meta-atom lattices with smaller filling factors.
Fig. S7 shows more examples of polaritonic energy spectra calculated for the same meta-atoms using smaller values of filling factor $f$. Expectedly, the Rabi splitting get smaller with decreasing filling factor as less and less meta-atoms occupy the same volume.
Fig. \ref{fig3b} shows the mode splitting $\Omega = 2g$ as a function of the center-to-center inter-particle distance $a/r$.
Even for the center-to-center inter-particle distance $a = 10 r$ corresponding to the moderate value of the filling factor $f = 0.0167$ the mode splitting in the ensemble of $r = 500$ nm meta-atoms reaches a sizeable fraction of the resonant energy, $\Omega_R \approx 0.2 \omega_{pl}$, and remains well above the characteristic non-radiative plasmon decay rate of 50 meV.



\subsection{Polaritonic spectra with spheroidal meta-atoms}

Next we analyze the behavior of the collective coupling constant in arrays of closely packed spheroidal meta-atoms, \ref{fig4}. This geometry is a good analytical approximation for elongated nanorods or nanodisks, which have been employed in a number of works studying strong and ultra-strong coupling in systems of meta-atoms \cite{bisht2018collective,Baranov2020}. Elongating one of the axes of the spherical meta-atoms affects their resonant properties and thus the collective coupling constant.
The filling factor (Eq. \ref{Eq_1_f}) of the lattice remains the same due to proportional scaling of the dimensions of the array. 
Violated spherical symmetry of a spheroidal meta-atom, however, couples orthogonal vector spherical harmonics of the electromagnetic field, which is why the quasi-normal modes of a spheroidal meta-atom cannot be determined from a single characteristic equation, as in Eq. \ref{Eq_3}. Instead, an infinite chain of coupled equations must be used \cite{asano1975light}, which presents an extremely excruciating problem.


For this reason, we obtain the eigenfrequency spectra of spheroidal meta-atoms numerically with the use of finite-element software COMSOL Multiphysics along with the specialized MAN package (Modal Analysis of Nanoresonators) \cite{wu2022modal}. The usage of additional software is desired because the in-built COMSOL eigenfrequency solver demonstrates insufficient convergence in problems with highly dispersive materials. We employ the QNMEig solver which implements auxiliary-field technique for finding  quasinormal modes (QNMs) \cite{yan2018rigorous}




Figures S8 and S9 show the resulting transition dipole moments and normalized collective coupling strengths of prolate spheroids as a function of longer semi-axis $a$ for a fixed plasma frequency $\omega_P = 10$ eV and varying aspect ratio (AR) (the special case of spherical meta-atoms is included as $AR = 1$). Overall, they demonstrate the behavior qualitatively similar to the ones obtained for spherical meta-atoms.
Similarly to spheres, normalized coupling constants quickly approaches an upper limit $\bar{\eta}$ with increasing longer semi-axis.
This upper limit, however, is evidently different for every aspect ratio. Given the behavior of this constant for spherical meta-atoms, we expect it to be universal for any value of plasma frequency and only depend on the geometric AR of meta-atoms.


To that end, we study how the upper limit of the normalized coupling strength $\bar{\eta}$ depends on the spheroid aspect ratio.
We approximately determine this limit as the maximal value of $\eta$ for each AR from our data points, which is a good measure given how quickly the normalized coupling strength approaches the plateau.
The resulting behavior of $\bar{\eta}$ shown in Fig. \ref{fig4} reveals a nearly linear dependence of the upper limit of normalized coupling with the meta-atom aspect ratio, in particular reproducing the ultimate value of $\bar{\eta} \approx 1.2$ obtained for spherical meta-atoms. This suggests that prolate plasmonic meta-atoms are more efficient for achieving deep strong coupling between light and matter.

The above analysis clearly suggests that prolate spheroids are beneficial for reaching higher values of normalized coupling strength $\eta$. However, elongating a metallic nanoparticle comes at a price of red-shifting its resonances.
Therefore, a natural question arises: what is the optimal shape (aspect ratio) of the metallic meta-atom that maximizes the absolute coupling constant $g_{\mathbf{k} }$ for a given resonant frequency $\omega_0$?

To that end, we utilize the same data and present the absolute value of the collective coupling constant in eV for all studied meta-atom aspect ratios (including the spherical case with $\textrm{AR} = 1$) as a function of the meta-atom resonant frequency, Fig. \ref{fig5}. The data allows to conclude that not only prolate spheroids offer high normalized coupling constant, but also yield the highest absolute value of the collective coupling constant for any given resonant energy.
This observation suggests that highly elongated metallic meta-atoms are, perhaps, the optimal geometry for the purpose of realizing collective polaritonic states with the largest coupling constant and Rabi splitting \cite{hertzog2021enhancing}.

\begin{figure}[t!]
\centering\includegraphics[width=1\columnwidth]{Fig_8.pdf}
\caption{Absolute values of the collective coupling strength in closely packed arrays of metallic spheroids as a function of the resonant meta-atom frequency $\omega_0$. The data is obtained for a fixed $\omega_p = 10$ eV.}
\label{fig5}
\end{figure}


Before concluding, we would like to emphasize the value of our results despite the number of crude simplifications and assumptions we have made in the model. Indeed, we assumed only electric dipole field-matter coupling, and neglected dipole-dipole inter-particle interactions in the total Hamiltonian. At the same time, all the values of the collective coupling strength have been obtained in the limit of densely packed particles, where the near-field inter-particle interactions may become crucial. 
However, the developed formalism easily allows one to recalculate the coupling strength for less dense arrays with $f \ll 1$, where the inter-particle interactions may be ignored.
In this case, all the quantities will be equally scaled by a factor of $\sqrt{f}$ (see Eq. \ref{Eq_9} and \ref{Eq_17}).
Therefore, the conclusion about the prolate meta-atoms yielding the highest coupling strength will remain true for any given filling factor of the meta-atom array, which is one the central findings of our work.


\section{Conclusion}
To conclude, we have studied collective polaritonic states formed by arranged plasmonic meta-atoms interacting with the free space optical field.
Almost linear scaling of the transition dipole moment of spherical meta-atom with its radius causes the collective coupling constant to quickly enter the ultrastrong and deep strong coupling regime before approaching a universal upper bound.
The resulting bound is universal for all plasma frequencies and is determined only by the geometry of the meta-atom.
The corresponding polaritonic energy spectra, calculated with the use of the developed Hamiltonian model, exhibit large values of Rabi splitting and polariton gaps.
Similar analysis of the arrays of spheroid meta-atoms showed that the normalized collective coupling constant and it upper bound increases with the aspect ratio of elongated metallic meta-atoms.
Furthermore, for any given resonant energy highly elongated spheroidal meta-atoms exhibit the highest absolute coupling constant.
These results should open up new prospects for realizing polariton states with artificial meta-atoms.
\\




\section{Acknowledgements}
Authors acknowledge fruitful discussion with Timur Shegai and Andrey Bogdanov. B.R. acknowledges fruitful discussions with Prof. Hans-Rudolf Jauslin and Prof. G\'erard Colas des Francs.
The authors gratefully acknowledge the financial support from the Ministry
of Science and Higher Education of the Russian Federation (Agreement No. 075-15-2022-1150). The Laboratoire Interdisciplinaire Carnot de Bourgogne is a member of EIPHI Graduate School (contract ANR- 17-EURE-0002).
D.G.B. acknowledges support from Russian Science Foundation (grant No. 21-72-00051) and BASIS Foundation (grant No. 22-1-3-2-1).

