Abundance of cavity-free polaritonic states in resonant materials and nanostructures

Strong coupling between various kinds of material excitations and optical modes has recently shown potential to modify chemical reaction rates in both excited and ground states. The ground-state modification in chemical reaction rates has usually been reported by coupling a vibrational mode of an organic molecule to the vacuum field of an external optical cavity, such as a planar Fabry-P\'erot microcavity made of two metallic mirrors. However, using an external cavity to form polaritonic states might: (i) limit the scope of possible applications of such systems, and (ii) be unnecessary. Here we highlight the possibility of using optical modes sustained by materials themselves to self-couple to their own electronic or vibrational resonances. By tracing the roots of the corresponding dispersion relations in the complex frequency plane, we show that electronic and vibrational polaritons are natural eigenstates of bulk and nanostructured resonant materials that require no external cavity. Several concrete examples, such as a slab of excitonic material and a spherical water droplet in vacuum are shown to reach the regime of such cavity-free self-strong coupling. The abundance of cavity-free polaritons in simple and natural structures questions their relevance and potential practical importance for the emerging field of polaritonic chemistry, exciton transport, and modified material properties.


I. INTRODUCTION
Strong coupling is a distinct regime of light-matter interactions, which is reached when resonant optical modes and material excitations (electronic, vibrational, etc.) exchange energy faster than they lose it to the environment. This fast energy exchange gives rise to new quasiparticles; polaritons. 1,2 Following the achievement of strong coupling, fundamental questions have arisen, such as whether polaritonic states could influence material and chemical properties. 3,4 The possibility of affecting chemical reactions in the excited state (i.e. photo-chemistry) seems to be quite well understood today. Several experimental studies show that strong coupling can affect photo-chemical processes [5][6][7] and these results are essentially in agreement with theoretical descriptions. [8][9][10][11][12][13][14] Thermally-activated ground-state chemical reactions in the vibrational strong coupling (VSC) regime seem to be more controversial. In this case, just a few experimental reports exist [15][16][17][18] and, there is far less conclusive agreement with theory. [19][20][21][22][23][24][25][26] An example of this controversy is given by common intuition, by stating that ground-state chemical reactions are local, i.e. the chemistry at position A does not depend on any far-away position B. Nevertheless, recent experimental a) Electronic mail: timurs@chalmers.se. studies [15][16][17][18] suggest that such a local approximation to the chemistry might not always hold. Contemporary empirical observations, 3,4,27 indicate that the chemistry can be affected by resonant, collective and non-local variables. For example, the reaction probability at a specific position depends on the far away presence (or absence) of metallic mirrors forming a Fabry-Pérot (FP) cavity. Specifically, it depends on the set of resonant optical modes that the mirrors form. This set of resonant optical modes and its relevance to chemistry is especially important, since they distinguish polaritonic chemistry from other common non-local but non-resonant contributions, such as screening, electrostatic, or solvent polarity effects. [28][29][30] Moreover, while standard chemical theory states that the reaction rate depends on the concentration of molecules, concentration-dependent collective Rabi splitting influences the rate further by modifying the energy levels. Although this train of thought is heavily debated in the community, mainly due to a lack of theoretical understanding of the recent experimental observations, we believe the potential of this new research frontier is quite high. Of course, further experiments and theoretical developments are needed to clarify this and other subtle issues.
Without trying to resolve the discrepancy of the groundstate chemistry in the VSC regime at this point, here we wish to focus on another potentially relevant issue. The so far experimentally reported observations of chemistry modifications under VSC exclusively use FP cavities comprised of two metallic mirrors. [15][16][17][18] Then, a natural question is if the role of the mirrors in these experiments is only to confine elec-tromagnetic modes in order to form polaritons, or if it goes beyond that. Thus the importance of discussing the concept of "cavity" in the strong coupling regime and understanding its relevance to polaritonic chemistry.
Traditionally, one requires an external cavity to reach strong coupling, such as a FP cavity, a photonic crystal cavity, or a plasmonic nanocavity, schematically visualized in Fig. 1(a). Assuming that the role of the cavity is solely to confine modes, we argue that using external cavities may limit the scope of potential applications, as such cavities are impractical. 31 There also exist a number of so-called "open" cavities, which might be more accessible and useful, 1,2,6 but also require an external object to provide the optical mode confinement. Here, we show that to achieve strong coupling, the cavity does not need to be external. Instead, electronic or vibrational excitations in bulk or nanostructured materials can self-couple to an optical mode sustained by its own geometry, as illustrated in Fig. 1(b). We consider a generic medium with a Lorentz dielectric function, motivated by the fact that optical properties of most materials can be modelled as such.
It is well known that objects with refractive index different from their surrounding can sustain localized optical eigenmodes. 32 The characteristics of the modes depend on the geometry and refractive indices of the involved objects and its surroundings. Once such an optical resonance is found at a proximity of the material resonance, these optical and material resonances can hybridize and, provided the oscillator strength of the transition is high, give rise to polaritons. The Rabi splitting in such a self-hybridized polariton can approach (but cannot exceed) the so-called bulk polariton splitting, which is a function of oscillator strength only and does not depend on parameters of the cavity. Such self-coupled or cavity-free polaritons have been realized in many relevant systems, such as slabs of hexagonal boron nitride (hBN), transition metal dichalcogenides, perovskites, and others. [33][34][35][36][37][38][39] However, their relevance for the polaritonic chemistry and modified material properties has not been raised in the community.
Originally, the problem of material resonance from the quantum optical light-matter coupling perspective has been solved by Hopfield in his seminal work in 1958. 40 Hopfield considered a three-dimensional (3D) continuous semiconductor and modelled it using a quantum coupled harmonic oscillator model, but did not consider materials shaped into a specific object(s) that can sustain a set of discrete well-defined electromagnetic modes.
As we show here, cavity-free polaritons are found in a great variety of structures, ranging from bulk materials to spherical micro droplets, Fig. 1(b). The sole existence and apparent abundance of such cavity-free polaritons questions their applicability in polaritonic chemistry and, more generally, in any polariton-induced modifications of materialrelated properties. 30

II. REVISITING STANDARD CAVITY-EMITTER POLARITONS
We start by revisiting the standard light-matter interaction Hamiltonian in Hopfield's formulation. 40 This formulation is essentially a coupled harmonic oscillator problem approached from the quantum optical perspective. We note that this is not exactly the same model as the famous Jaynes-Cummings or quantum Rabi models, which deal with quantum emitters (QEs) coupled to a cavity mode. However, within the weak excitation regime, this formulation is sufficiently appropriate even for QEs, and even more within the scope of this study since we focus on macroscopic Lorentz-resonant descriptions.
The full Hopfield Hamiltonian for a cavity mode oscillating at a frequency ω c interacting at a rate g with a material transition at a frequency ω 0 , written in the Coulomb gauge and including the diamagnetic term, reads: whereb andâ are the annihilation operators of the material transition and the cavity mode, respectively, and g is the coupling constant. 41 Two eigenvalues (resonant transitions) of this Hamiltonian are given by: With the material oscillators homogeneously occupying the cavity mode volume V , the coupling strength g = µ ρV /3E vac ω 0 /ω c (with the factor 1/3 accounting for their isotropic orientation). Thus, the coupling strength depends on the transition dipole moment of the material transition, µ, and the volume density of these dipoles (oscillators), ρ as well as the vacuum electric field of the cavity, E vac = hω c /2ε ∞ ε 0 V . When the coupling strength is large enough, two polaritons with different energies (Eq. 2) are formed, which has been explained in great detail before. 2,42 From the above description, one can deduce that two main ingredients are necessary for achieving strong coupling: an optical cavity with a large quality factor (Q) and a small mode volume (V ) and an electronic (or vibrational) transition with a high transition dipole moment. Following such suggestions, there has been extensive research to obtain the best optical cavities, ranging from high-finesse FP cavities composed of distributed Bragg reflector (DBR) mirrors to plasmonic nanocavities. Many of these cavities are covered by Baranov et al. 2 and some are depicted in Figure 1.
In the examples above, we see that a stand-alone cavity is important for realizing the standard cavity-quantum electrodynamics (QED) scenario, that is, the case of a single-QE strongly coupled to an external optical cavity. Such a scenario is challenging to realize with the cavity-free approach, described below. Moreover, a single QE problem can not be treated in terms of a macroscopic Lorentz dielectric functions even with the Hopfield Hamiltonian, as we do below, but instead requires a true cavity-QED approach, which we do not focus on here (these approaches are extensively covered in the literature 1,43 ).

III. SELF-HYBRIDIZED POLARITONS
An optical mode is a concept that is not limited to the optical resonators and cavities discussed above. An optical mode is a solution of the source-less Maxwell equations in a given geometry. Thus, a cavity can be made of any material, including the material that makes up the relevant electronic or vibrational transitions themselves. These are precisely the scenarios we focus on in the following sections.
A. 3D case: bulk polaritons Let us start by revisiting bulk polaritons, hybrid lightmatter states of an unbounded dielectric medium homogeneously filled with resonant transitions. This situation has been described by Hopfield in the celebrated 1958 paper 40 (see also Mills et al. 44 ). Instead of a localized eigenmode of a cavity, the photonic modes in this scenario are represented by plane waves propagating in the unbounded dielectric medium with permittivity ε ∞ . Those waves have a continuous spectrum ω = kc/n with k being the vacuum wave vector, c the speed of light, and n = √ ε ∞ the refractive index of the medium.
In the classical electromagnetic formalism, an isotropic resonant medium can be described by a Lorentzian permittivity: with ω P = ρe 2 /3ε 0 m being the plasma frequency, where ρ is the volume density of oscillators and 1/3 accounts for their isotropic orientation, e and m are, respectively, the electron charge and mass and ω 0 and γ are the resonance frequency and linewidth, and f the oscillator strength. The latter is expressed via the microscopic parameters of the medium as f = 2 mω 0 e 2h |µ| 2 , where µ is the transition dipole moment of the material oscillator (see Methods). For simplicity, in Eq.
(3) and throughout this study (except for the case of hBN), we deal with only one resonance with an isotropic random orientation of transitions. More complicated cases, including multiple Lorentz resonances and anisotropy, naturally follow from this simplified model.
Bulk polaritons correspond to roots of the dispersion equation: One can look for solutions of this equation either with complex ω and real k, or complex k and real ω. These two types are both appropriate solutions of Maxwell's equations, but suited for different purposes. 45 We are interested in complexω solutions of the dispersion equation, since they better reflect the transient decay of hybrid light-matter states of the coupled system.
Quantitative characteristics of the problem, such as the vacuum Rabi splitting and the magnitude of the polaritonic gap, can be obtained from the Hamiltonian formulation of the problem that neglects the exciton dissipation. 40 Following the standard procedure of quantizing the electromagnetic field in free space, let us consider a quantization box of volume L 3 with periodic boundary conditions. 46 The vacuum electric field of a photon with frequency ω c reads E vac = hω c /2ε ∞ ε 0 L 3 . Recalling the expression for the oscillator strength f and plasma frequency ω P of a Lorentz medium, we express the resulting coupling strength as g C = (ω P /2) f ω 0 /ε ∞ ω c . The vacuum Rabi splitting (with zero losses) is then obtained as the difference between the two polariton energies (Eq. 2) for the zero-detuned photon (ω c = ω 0 ): where g 0 = (ω P /2) f /ε ∞ is the zero-detuning coupling strength.
The upper edge of the polariton gap is obtained by calculating the upper polariton energy in the limit k = 0 and is ω 0 = ω 2 0 + 4g 2 0 , whereω 0 is the re-normalized resonance frequency. The lower edge of the gap is found by calculating the lower polariton energy in the limit k → ∞, and is exactly the uncoupled exciton energy ω 0 . This results in the polaritonic gap which simplifies to ∆ pol ≈ 2g 2 0 /ω 0 under the assumption of 4g 2 0 /ω 2 0 ≪ 1, that holds in the standard strong coupling regime.
The presence of losses in the dielectric formalism makes the eigenfrequencies complex and changes their dispersion, however, the original Hopfield model of bulk polaritons does not include any dissipation channels. Subsequent works extended Hopfield's model to account for exciton decay by considering coupling between them and a continuum bath of harmonic oscillators. 47 Such a model predicts that bulk polaritons appear if and only if the condition holds. It is easy to see that this condition is equivalent to (ω P /2) f /ε ∞ = g 0 > γ/4, which is the strong coupling criterion in the standard Jaynes-Cummings model with a lossless cavity mode. 47,48 Thus, by knowing the parameters of the material resonance in Eq. (3), it is possible to estimate the bulk coupling strength at zero detuning g 0 and compare it to γ/4 in order to conclude whether or not this specific material supports bulk polaritons. Figure 2(a-c) shows the dispersion of bulk polaritons (k vs real part of ω) obtained for a medium described by the permittivity in Eq. (3), with typical values for TDBC J-aggregates (ω 0 = 2.11 eV, γ = 0.1 eV, ε ∞ = 2.15) 49 , and varying oscillator strengths. First, for f ω 2 P = 0.0045 eV 2 , corresponding to a strongly diluted J-aggregate, we observe the mode attraction typical for the weak coupling regime, Fig. 2(a). In this regime the energies of the eigenmodes of the medium slightly deviate from the bare values, but their dispersions cross near the zero detuning region.
For a larger oscillator strength f ω 2 P = 0.445 eV 2 , see e.g. Balasubrahmaniyam et al., 49 one can clearly see an anticrossing between the two polaritonic modes, with the vacuum Rabi splitting around 455 meV. We calculate the vacuum Rabi splitting as the energy separation along the real frequency axis between polaritonic modes obtained at a zero detuning between the electronic resonance of the material and the bare cavity mode (in turn, obtained by setting f = 0). This choice corresponds to nearly 50%/50% excitonic/cavity polariton wave-function amplitudes, and we use this definition throughout the text. On the same plot we also show dispersion of complex-k/real-ω solutions (in a solid red line), found simply as k = ε(ω)ω/c. The two solutions coincide away from the resonance, but differ in the vicinity of ω 0 . While the two branches of complex-ω solutions are disconnected, the two branches of the complex-k solutions are connected by the region of anomalous dispersion.
For even larger oscillator strength, f ω 2 P = 3.116 eV 2 (we note here that this value is unrealistically large for Jaggregates and we give it here only for the sake of theoretical argument; there may be other material systems, where such oscillator strengths may be reached), the dispersion shows a larger vacuum Rabi splitting of about 1.2 eV and also displays a polaritonic gap -a region of frequencies with no propagating modes within it, shown in Fig. 2(c). 40,43,50 The corresponding complex-k solutions, at the same time, span over the whole range of frequencies, but become highly damped within the polariton gap by acquiring a large imaginary part of the wave vector k (see Fig. S1).
In addition to the anticrossing in the k − ω plane, it is instructive to inspect trajectories of the roots of the dispersion in Eq. (4) in the complex-ω plane. Fig. 2(d-f) show trajectories of the complex frequency roots for the three cases of weak, strong, and ultrastrong coupling. In the weak coupling, f ω 2 P < ε ∞ (γ/2) 2 , the "photonic" root of the equation crosses the exciton position and acquires a small imaginary component in the zero detuning region, while the "excitonic" root (the one originating at the complex exciton frequency, marked with a yellow star in Fig.2) moves along a finite curve around ω 0 − iγ/2. However, in the strong and ultrastrong coupling regimes the low-energy solution terminates exactly at the exciton energy, while the upper solution originates on the other side of the polaritonic gap, Fig. 2(e,f).
In closing this part we want to emphasize that bulk polaritonic states and the whole anticrossing picture exist in an optically large piece of a resonant material without any external cavity. The examples given above clearly illustrate that a large (compared to the resonant wavelength) unstructured excitonic or phononic material already possesses polaritonic modes, provided that the oscillator strength of the material exceeds the threshold value. Reducing the dimension of the system (by breaking the translational invariance in one, two, or all three dimensions) does not break this picture: the reduced system still possesses optical modes (either complex-frequency resonant states, or real-frequency guided modes) that couple with the material resonance.
Let us first reduce the bulk material in one dimension (along the z-axis, to be specific) and consider a slab with a Lorentzian permittivity ε(ω) in vacuum. The eigenmodes of this system can be quantified by the (real) in-plane momentum k x , k y . 51 In the absence of the material resonance, the modes above the light-line are FP resonances with complex frequencies (not to be confused with real-frequency leaky modes); 52 the modes below the light-line are proper waveguide modes with real frequencies. 53 Upon inclusion of the material resonance, it couples to modes both below and above the lightline. Rabi splitting previously mentioned.
TE-polarized eigenmodes of a 200 nm thick J-aggregate slab with a higher oscillator strength f ω 2 P = 3.16 eV 2 demonstrates anticrossing both below and above the light-line, Fig.  3(b). This J-aggregate was analyzed in Fig. 2(c,f), and corresponds to the bulk ultrastrong coupling case. The location of the Rabi splitting with TE 1 is marked with a purple arrow, but its magnitude cannot be obtained because the upper polariton has a cut-off given by the light-line and the corresponding guided eigenmode does not exist.
Next, we demonstrate self-hybridized polaritons in another practically relevant system: a slab of hexagonal boron nitride (hBN). hBN exhibits two prominent vibrational modes at 169 and 95 meV, respectively. 54 We will focus on TEpolarized modes, which couple only with the in-plane vibrations of hBN at 169 meV. The in-plane permittivity component of hBN can be described by a single Lorentz term with ε ∞ = 4.45, f ω 2 P = 0.064 eV 2 , and γ = 1 meV, which approximates the experimental data in the relevant spectral range ( Fig.  S1(g)). 54 Fig. 3(c) shows eigenfrequencies of TE-polarized eigenmodes of a 1.75 µm thick hBN slab. Clearly, there are anticrossings both below (• and ▽) and above ( and △) the light-line thanks to the large oscillator strength of hBN's vi-brational transition. For the same reason, we also observe the lower polariton branch of the 2 nd FP mode above the lightline, shown in blue △, and lower polariton of the TE 2 waveguide mode below the light-line (orange •). The Rabi splitting of the TE 1 waveguide mode, which occurs at about k x = 1.5 µm −1 , reaches a value of ≈ 116 meV, which is also below the bulk Rabi splitting, 2g 0 ≈ 120 meV. Remarkably, below the light-line, the two fundamental modes -TE 1 and TE 2 -show pronounced anticrossing, even though the TE 2 mode dispersion in a dielectric slab without the exciton does not approach the exciton energy (see Fig. 3(c)).
This behavior becomes more evident when we look at trajectories of the eigenmodes' poles of the slab in the complexω plane. It turns out, that even for a vanishingly small oscillator strength f ω 2 P there are countably many poles in an arbitrarily small vicinity of the complex exciton frequency ω 0 − iγ/2 (see Fig. S3(b)), where the Lorentz permittivity ε(ω) has a simple pole (shown in Fig. S3(c)). 55,56 Each of these poles represents an eigenmode of the Lorentzian slab: either a radiative mode, or a localized waveguide mode, depending on the position of the pole frequency with respect to the complex plane cut.
The origin of this cluster of poles can be appreciated by considering the Lorentzian permittivity ε at complex frequencies, Fig. S3(c). One can clearly see that just below the com-plex exciton frequency, ω 0 − iγ/2, the permittivity takes arbitrarily large values, thus, in theory, allowing for the existence of higher-order FP and waveguide modes, that are otherwise inaccessible at real frequencies, where the permittivity is bounded by |ε(ω)| ≤ ε ∞ + f ω 2 P /γω 0 . Because these poles are located so close, they usually merge into a single spectral feature in real-frequency spectral responses, often attributed to "uncoupled molecules". 57,58 In fact, as the complex frequency plane reveals, this feature stems from a whole new multitude of eigenmodes located in the small vicinity of the complex exciton energy (Fig. S3). Of course, only a finite number of these eigenmodes are physical, since the permittivity becomes highly non-local and loses its physical meaning as soon as the wavelength inside the material λ / Re( ε(ω)) at complex ω approaches the microscopic length scale of the medium near the pole.
Complex-plane trajectories also clearly reveal the Rabi splitting. Above the light-line trajectories of the eigenmodes are limited by the low frequency cut-off on the one side, and by the light-line on the other side, and for this reason are difficult to interpret (Fig. S4). Below the light-line, however, only the upper polariton branches are bounded by the light-line on the low energy side, Fig. 3(d), and the resulting trajectories are very similar to those for bulk polaritons (see Fig. 2(e)). Qualitatively similar pictures of anticrossing in the ω − k x space and in the complex-ω plane are found for TM-polarized eigenmodes (Fig. S5).
Since the waveguide modes of a dielectric slab have virtually infinite lifetime in contrast to radiative FP modes, anticrossing below the light-line appears at a much smaller oscillator strength. For an electrically thin film we can establish an approximate analytical criterion of strong coupling for waveguide modes. Assuming an electrically thin dielectric film without the Lorentzian transition (k 0 L √ ε ∞ ≪ 1, where L is the thickness of the film), the dispersion of the TE 1 waveguide is approximately with δ k being a small parameter for a thin film (see Methods). The vacuum field of the bare TE waveguide mode in the center of the dielectric slab can be estimated as is the waveguide mode's effective width. Multiplying the vacuum field by the square root of the number of oscillators ρa 2 L in the volume a 2 L and the dipole moment µ, we find the coupling strength (at zero detuning, ω c = ω 0 , k 0 = ω 0 /c): Comparing the above g to γ/4, we find the threshold oscillator strength required for anticrossing below the light-line in a dielectric film: This condition may help estimate the minimal thickness L of an excitonic film required for the emergence of polaritonic modes below the light-line for a given oscillator strength.
C. 1D case: polaritons in long circular cylinders By reducing the system in one more dimension (say, along the y-axis) we arrive at one possessing translational invariance only in a single direction (along the x-axis). Thanks to this symmetry, eigenmodes of the one-dimensional system are characterized by a real-valued momentum k x . Similarly to the 2D case, eigenmodes with k x > ω/c are guided modes with real frequencies; eigenmodes with k x < ω/c are radiating resonances with complex frequencies.
For simplicity, let us focus on 1D systems with a circular cross-section. The states with k x = 0 (normal-incidence resonances of the cylinder) are classified by their polarization state (TE or TM) and the azimuthal number m. 59 Eigenmodes with k x = 0 split into TE and TM-polarized solutions only for m = 0 (the monopole harmonic). 60 For m > 0, TE and TM polarizations couple, giving rise to hybrid modes notated EHnm and HEnm. The HE11 mode (the dipole mode) is the only mode of a circular cylindrical dielectric waveguide not having a cut-off 60 (see ▽ in Fig. S6(a,c)).
As a practical example, we demonstrate self-hybridized polaritons in infinite circular cylinders of a perovskite CsPbCl 3 . This material is described by a Lorentz permittivity with ε ∞ = 3.7, ω 0 = 3 eV, γ = 87 meV, and f ω 2 P = 0.54 eV 2 that approximates the experimental data 36 (see Fig. S1(a)) . Polaritons and polaritonic emission in long whiskers made of perovskites of various families have been demonstrated in a number of works. 35,61 Fig. 4 presents frequencies of TM-polarized monopole (m = 0) eigenmodes for CsPbCl 3 infinitely long circular cylinder with a radius of 250 nm. The resulting spectrum clearly shows anticrossing above the light-line (involving the 2 nd order radiating mode) and below the light-line (involving the TM 0,1 waveguide mode). The anticrossing below the lightline features a Rabi splitting of about 312 meV, which compares to bulk Rabi splitting of CsPbCl 3 of 382 meV. Complexω trajectories of eigenfrequencies exhibit a familiar behavior with upward frequency pulling to the real axis above the lightline, and downward pulling below the light-line, Fig. 4(b,c). Qualitatively similar eigenfrequencies are obtained for TEpolarized monopole modes (Fig. S6), and for the dipole TEM modes (Fig. S7).

D. 0D case: polaritons in spheres
Finally, if we reduce the system in one more dimension, we end up with a compact scatterer supporting a discreet spec- trum of localized eigenmodes. 62,63 Because of the lack of translational invariance in any direction, all eigenmodes of the scatterer are complex-frequency radiating solutions. 64 To keep the analysis simple, we consider the most symmetric compact scatterer: a sphere, whose eigenmodes are the well-known Mie resonances. 65 These resonances split into TE-and TM-polarized modes, which are further classified by the angular index l, and the radial number N. l = 1 modes are usually referred to as the dipole ones, l = 2 as the quadrupole ones, and so on. A detailed description of whispering gallery exciton-polaritons in GaAs spheres has been presented by Platts et al. 33 Because of the lack of translational invariance, the anticrossing in spheres cannot be seen by scanning the eigenmode's momentum k. Instead, one usually varies the radius of the sphere, what causes the eigenmodes to move in the ω-space.
As a practically interesting example of self-hybridized polaritons in spheres, we consider liquid water droplets. Liquid water has a strong vibrational resonance around 420 meV with the oscillator strength of f ω 2 P = 0.035 eV 2 . 66 Therefore, one may naturally expect formation of vibrational polaritons in Mie resonant micron-sized water droplets. We describe water by a Lorentz permittivity with ε ∞ = 1.75, ω 0 = 0.42 eV, γ = 0.048 eV, and f ω 2 P = 0.035 eV 2 that approximates the experimental data, see Fig. S1(d). 66 Fig. 5(a) shows eigenfrequencies of TE 1,1 (magnetic dipole) eigenmodes of a water sphere as function of its radius r. The trajectory of the TE 1,1 eigenmode (the first index indicates the angular number of the mode l, and the second index stands for the radial number of the mode N) shows a clear anticrossing with a magnitude of about 30 meV for a sphere with r ≈ 1.1 µm (Fig. 5(c)). This value is significantly lower than the corresponding bulk Rabi splitting for water ω P f /ε ∞ ≈ 100 meV due to high radiative loss of the TE 1,1 mode. The small value of the Rabi splitting in this case is also manifested in the shape of the dispersion curves, similar to the exceptional point (EP) regime with the square root dependence of eigenenergies on the perturbation parameter (such as the radius). 67 We find that 1.1 µm is the smallest radius for which polaritonic modes appear in water spheres. Although the TM 1,1 (electric dipole) eigenmode of smaller r ≈ 600 nm spheres with ε ∞ = 1.75 matches the exciton position, it has much lower Q-factor, which prevents strong coupling in spheres of this size (Fig. S8).
Nevertheless, the TM 1,2 mode of r ≈ 1.7 µm spheres does demonstrate a clear anticrossing with the water vibration, Fig.  5(b,d), producing a Rabi splitting of about 70 meV, which approaches the bulk splitting of 100 meV. Complex-ω trajectories in both cases reveal the characteristic eigenfrequency pulling to the real axis, Fig. 5(c,d), originating from large radiative losses of the bare optical mode compared to the exciton linewidth. With increasing radius of the sphere new eigenmodes corresponding to higher orbital and radial numbers (and having higher Q-factor) will match the vibrational resonance position and show additional anticrossings and polaritonic modes. We further investigate how the attainable Rabi splitting in water spheres depends on the eigenmode polarization and the orbital number. Fig. 6(a) shows the Rabi splitting for TE and TM eigenmodes as a function of the orbital number l. One can see that the Rabi splitting of TE modes systematically exceeds that of TM modes for all orbital numbers, but both gradually approach the bulk water splitting, 2g 0 , which for water is ∼ 100 meV (see table S1).
Coincidentally, micron-sized droplets of liquid water that support vibrational self-polaritons are quite ubiquitous in nature. For example, a common size of water droplets in fogs, mists, clouds, or steams falls exactly into the range of 1-10 µm that is needed for the formation of vibrational Mie selfpolaritons. Due to the surface tension these droplets adopt spherical geometry, thus forming an ideal natural platform for observation of vibrational self-polaritons in strong and even ultrastrong coupling regime. Note, that the bulk polariton splitting of liquid water falls into the vibrational ultrastrong coupling regime, since g H 2 O 0 ≈ 0.05 eV and η H 2 O 0 ≈ 0.12 > 0.1, which satisfies the ultrastrong coupling criterion. 26

A. Critical size for the cavity-free polariton formation
In this section, we turn our attention to the apparent fact that nanostructures supporting cavity-free polaritons have a certain critical size, below which no polaritonic behavior is observed. For example, water spheres below ≈ 1 µm in radius will not support polaritons. Such size-limiting behavior is, of course, not universal to water and should be observed in all Lorentz materials (e.g. r min ≈ 98 nm for perovskite and ≈ 180 nm for J-aggregates in Fig.6(a)). Neither is it exclusive to the spherical geometry. For dielectric slabs, where the fundamental waveguide modes do not have a cut-off, the physical reason for the critical size is the delocalization of the electric field of the mode and a resulting weak coupling strength. For circular cylinders (and their monopole modes) and spheres the size bound stems from the fact that the fundamental eigenmodes of electrically small systems are far blue-detuned with respect to the material transition resonance, in which case efficient hybridization is not possible. Finally, we would like to mention, that in the opposite case of exceedingly large 0D-2D objects, the result is trivial and reduces to the bulk case, which has no critical upper limit.
From the material-science perspective the existence of the critical size suggests that the polaritonic behavior and hence modified material properties are supported only within a certain size range, given by the material resonant properties, its surrounding(s), and its shape. From the computational point of view, it is also interesting to note that the particles of about the critical size, although are on the order of several tens to several hundreds of nanometers, contain a large number of oscillators (atoms or molecules). Thus, it is challenging to selfconsistently capture these collective cavity-free polaritons using atomistic methods, such as density functional theory 23,68 , due to their high computational cost.

B. The limit of Rabi splitting
The diversity of realizations of strong coupling, both using cavities (Fig. 1) and in the cavity-free format (Figs. 2-5), raises a question about the limit of g and Ω R that one can observe using Lorentz materials. A particularly important question is whether it is the material or the cavity, which plays the decisive role in this limit.
To shed light on this issue, we compare the bulk polariton coupling strength to the other realizations for the specific Lorentz material considered above, Fig. 6. From the presented data it appears that the lossless bulk polariton splitting, 2g 0 = ω P f /ε ∞ (see table S1 for the values for each material), indeed establishes an ultimate limit on vacuum Rabi splitting achievable with a given Lorentz material with any cavity mode. Indeed, a cavity of mode volume V supporting a localized eigenmode can host only as many as N = V ρ transitions. Assuming they all couple to the mode with an identical coupling strength, which is certainly an overestimation, the resulting Rabi splitting is seen to be bound from top by the bulk polariton coupling. The lossless bulk polariton splitting is in turn limited from above by setting f =ε ∞ =1, which gives Ω max R = 2g max 0 = ω P . Thus, the maximum bulk polariton splitting is principally limited from above by the volume density of oscillatorsρ (and a combination of natural constants).
It is also instructive to consider the maximum ratio between the coupling strength and the resonance frequency η max = g max 0 /ω 0 = ω P /2ω 0 . This ratio plays an important role, when the coupling strength approaches the transition frequency. Specifically, when η > 0.1, one usually distinguishes the ultrastrong coupling (USC) regime, while when η > 1 − the deep strong coupling (DSC) regime. Taking in mind the principal limits of the coupling strength, the parameter η is seen to be intrinsically related to the ratio of plasma and resonant frequencies in Lorentz materials. This ratio can be both below one and significantly above one (in principle unbound), depending on the material. However, one should keep in mind that in the limit of DSC the light and matter components of the polaritonic states actually decouple, presenting nearly uncoupled optical mode and material excitation as the eigenstates of the system. 69 This decoupling can be seen as the result of the screening by the diamagnetic term. In order to maintain equal cavity/matter fractions, one would have to tune the matter transition in resonance with the renormalized cavity mode. The same decoupling occurs in the dipole gauge, where the matter resonance experiences a depolarization shift due to dipole-dipole interaction of the electronic subsystem. 70 In fact, considering the ratio between the coupling constant and the renormalized resonance frequency η = g 0 /ω 0 = g 0 / ω 2 0 + 4g 2 0 , one can see that the latter is principally bound by η ≤ 1/2. 71 This hard limit indicates impossibility of reaching DSC with Lorentz media (if normalization is performed with respect to renormalized frequency).
Since the Rabi splitting is limited from top by the bulk polariton coupling strength, having a separate cavity at best allows confining polaritons in space, but does not increase the coupling strength due to the above limitation. As we showed above, the optical mode confinement can be attained without a separate cavity in the first place, with the exciton (of phonon) material itself playing the role of a confining resonator.
The bulk polariton splitting limit may serve as a practical and useful guidance for the selection of material platform in strong coupling experiments. This limit is the function of the material only, and does not depend on the nature of the cavity. For any practical realization of strong coupling the observed Rabi splitting can only approach that of the bulk, 2g 0 = ω P f /ε ∞ , an estimation of which is straightforward (requires no complicated simulation). Thus, a material with high enough oscillator strength will naturally support polaritonic eigenstates in all studied 0D-3D regimes.

V. CONCLUSION
To conclude, we have shown that polaritonic states are natural and ubiquitous to bulk materials and nanostructures that can be described by a generic Lorentz resonance(s). Such polaritons are universally encountered in systems that do not involve any stand-alone cavities or metallic mirrors. In these cavity-free systems, boundaries of the Lorentz medium play the role of the mirrors that allow formation of well-defined optical modes, which in turn couple to the resonant transition. We demonstrated such cavity-free polaritonic states in 3D-bulk materials, 2D-slabs, 1D-cylinders, and 0D-spheres. Concrete examples of real material systems, including spherical water droplets, perovskite nanocylinders, as well as slabs of J-aggregates and hBN were provided.
These findings might have important implications for the recently observed modification of material and chemical properties seen in electronic and vibrational strongly coupled systems. In particular, from the classical electromagnetism point of view, as well as the Hopfield Hamiltonian treatment, the cavity-free polaritons are no different from the more standard microcavity polaritons containing dense molecular layers. Thus, all observations relevant for the cavity polaritons could be just as well applicable to the cavity-free polaritons presented here, provided it is indeed the polaritonic nature of those modes, which is responsible for chemistry and material science modifications (other reasons, such as usage of highly polarizable metallic mirrors as was done in nearly all experimental studies so far might also play a role 31 ). However, in the cavity-free case, the exact polaritonic eigenmode picture depends on the size and shape of the resonant material and its environment. Thus, if we accept the hypothesis that polaritonic states themselves (and not the presence of metallic mirrors) are able to modify material and chemical properties, then these material and chemical properties should be dependent on the size and shape of this material, in agreement with the cavity-free polaritonic eigenmode picture. Clearly, such a cavity-free approach to modified material properties is beneficial from the practical point of view, as many of these structures exist naturally and are thus trivial to fabricate. A common concrete example considered here is the water droplets (encountered in fogs, mists, and clouds), which have a sharp cut-off below the radius of about 1 µm. Smaller water droplets will not support vibrational Mie self-polaritons and thus may experience different chemical properties in comparison to bigger droplets that support vibrational Mie self-polaritons and to bulk water that supports bulk vibrational polaritons. Similar consequences may apply to J-aggregates, hBN, etc., as well as to any Lorentz material with high enough oscillator strength of the relevant (electronic or vibrational) transition. This reasoning applies also to polariton-assisted modification of other material-related properties. We look forward to experimental tests of these predictions.

A. Permittivity of bulk medium
The permittivity of a resonant medium containing electronic or vibrational transitions can be obtained by calculating the electric dipole polarizability of a single such transition, and combining it with the volume density of homogeneously distributed transitions ρ in the medium.
The dipolar polarizability of a two-level system with the transition dipole moment µ in a weak external field can be written as : wheref is the transition's oscillator strength, m is the electron mass, e is the electron charge, and ⊗ is the outer product 72 (note ε 0 in the expression for polarizability, as required in SI units). Combining this with the volume density of the twolevel systems ρ and assuming random orientation of transition dipole moments, one easily obtains the expression for the polarization density and permittivity of the medium: where ε ∞ is the non-resonant background permittivity, ω 2 P = ρe 2 /(3ε 0 m) is the plasma frequency, and the factor 1/3 accounts for random isotropic orientation of dipoles in the medium.

B. Eigenfrequencies of a slab
Eigenfrequencies of a planar slab of thickness L are found as poles of the reflection coefficient in the complex frequency plane. The reflection coefficient of a TE or TM polarized plane wave with wavenumber k 0 = ω/c and in-plane momentum k x reads where with k z,1 = k 2 0 − k 2 x , k z,2 = ε(ω)k 2 0 − k 2 x being the zcomponents of the wave vector in vacuum and dielectric, respectively. k z,1 has a branching point at k x = k 0 , and a cut needs to be made. The field of the eigenmode outside the slab is described by the phase factor e ±ik z,1 z , where the plus and minus sign is applied in the regions z > L/2 and z < −L/2, correspondingly. The radiating modes above the light-line should have Re(k z,1 ) > 0, while proper waveguide modes have Im(k z,1 ) > 0. Therefore, we make the branch cut along the negative imaginary axis through ω = −i∞ and choose the Riemann sheet of the square root function such that Rek z,1 > 0 above the light-line, and Im(k z,1 ) > 0 below the light-line. This choice ensures that we find radiating eigenmodes and proper (localized) waveguide modes. Choosing another Riemann sheet of the square root in the definition of k z,1 results in dispersion of improper (diverging) modes below the lightline. The choice of sheet for k z,2 is not important, because the slab contains both waves with ±k z,2 .
C. Dispersion of the TE 1 waveguide mode in a thin dielectric film Dispersion of TE waveguide modes in a dielectric film is given by the roots of the characteristic equation 1 − r 2 TE e 2ik z,2 L = 0. This is a transcendental equation and in a general case cannot be solved analytically. However, in the case of an electrically thin film (k 0 L √ ε ∞ ≪ 1) the equation can be linearized and an approximate solution of the fundamental TE 1 mode can be found.
For an electrically thin film, dispersion of the TE 1 mode is sandwiched between the light-line of the environment (k 0 = ω/c) and that of the slab (k = √ ε ∞ ω/c), and follows very close to the environment light-line; therefore, we can write the propagation constant of the guided mode as k x = k 0 + δ k, where δ k ≪ k 0 . The z components of the wave vectors in the environment and in the slab can be expanded as follows: Note that k z,1 scales as O( √ δ k), whereas the leading correction to k z,2 scales as O(δ k). The exponential phase factor can be expanded in a Tailor series e 2ik z,2 L ≈ 1 + 2ik z,2 L. Plugging these expansions into the characteristic equation, we obtain: Keeping only the leading power O( √ δ k) of the k x expansion and the constant term O(1) in this equation, we find δ k = k 3 0 L 2 (ε − 1)/8, and the approximate dispersion relation of the TE 1 , correspondingly, can be written as: D. Eigenfrequencies of a cylinder TE (TM) modes of an infinitely long cylinder are defined as those having their electric (magnetic) field strictly perpendicular to the cylinder axis. Eigenfrequencies of TE01 and TM01 (monopole transverse electric and transverse magnetic) waveguide modes of an infinitely long circular cylinder of radius a are found as roots of the characteristic equation. 59,60 where J 1 (z) and H 1 (z) = H (1) 1 (z) are Bessel and Hankel functions of the first kind, respectively, k ρ,1 = k 2 0 − k 2 x and k ρ,2 = ε(ω)k 2 0 − k 2 x are the radial wavenumber of the eigenmode in vacuum and dielectric, respectively.
The definition of k ρ,1 requires a branch cut at k ρ,1 = k 0 . Radial dependence of the eigenmode field outside the waveguide is given by the factor H 1 (k ρ,1 ρ). Therefore, making the branch cut going along the negative imaginary axis through ω = −i∞ and choosing the Riemann sheet yielding Re(k ρ,1 ) > 0 above the light-line results in the spectrum of radiating eigenmodes and localized waveguide modes. 59 Eigenfrequencies of the HE11 (dipole) mode are found as  In this section we explain how we obtain the poles in the complex−ω plane, we show the fittings for the permittivities used throughout the text, the bulk polariton dispersion with the complex−k and the values of the bulk Rabi splitting, 2g 0 , for all the materials used.
First, the video in the supporting material shows the poles in the complex−ω plane of a bulk TDBC J-aggregate with a permittivity given by ω 0 = 2.11 eV, γ = 0.1 eV, ε ∞ = 2.15 and f ω 2 P = 0.445 eV 2 , ? for different values of k from 5 to 25 µm −1 . The position of the poles, for each k gives us the information of both dispersion -k vs Re(ω-and trajectories of the poles in the complex-ω plane -Re(ω) vs Im(ω) -of the eigenmodes shown in Fig. 2(b,e). The value of the intensity is not considered in the main text, but it relates to the excitation rate of the specific mode. The video shows a clear Rabi splitting at k = 16µm −1 and depicts how the lower and the upper polariton are separated by the polariton gap, which can be found as defined in equation 6 in the main text.
The first column of Fig. S1 shows, in dashed lines, the data of the permittivity of (a) perovskites CsPbCl 3 , (d) water and (g) hBN. As well as, the Lorentz model fitting for each material, in solid lines, which was used for the calculations. Note that all the frequencies ω in the labels correspond tohω since the units used throughout the paper are eV.
The second column of Fig. S1 shows the bulk polariton dispersion by considering complex−k (k = ε(ω)ω/c) for hBN, perovskites and water.The same information for J-aggregates is found in the upper panels of Fig. 2 in the main text. We can also see the two branches of each bulk polariton (b,e,h) connected by a region of anomalous dispersion. Finally, the last column shows the damping of the solutions within the polariton gap due to the large values of Im(k).

II. 2D: POLARITONS IN LORENTZ SLABS
This section shows the slabs calculations for the TE trajectories of the radiating modes, the uncoupled optical modes, poles in the complex−ω plane and the calculations for hybridization with TM modes. Figure S2 shows the dispersion of the uncoupled optical modes above (in blue) and below (in orange) the light-line (obtained by setting f = 0). The poles are plotted on top of the colormap for the reflectivity of the slab. The first two panels correspond to a 200 nm thick J-aggregate slab with (a) for TE polarization and (b) for TM. (c) shows the same information for the TE modes of a 1.75 µm hBN slab. Such information is necessary to find k x such that there is zero detuning between the resonance and the optical mode. Figure S3 shows the poles in the complex−ω plane (similar to the ones in the supporting video) for a 200 nm thick J-aggregate slab f ω 2 P = 3.116 eV 2 at normal incidence, k x = 0, thus all of the poles are radiating modes corresponding to the upper polariton (UP), the lower polariton (LP) and the 0 th -order FP mode below the exciton gap. The countable many poles in the vicinity of ω 0 − iγ/2, are described carefully in the main text. In short, they correspond to higher order modes given by the steep increase of ε(ω), due to the pole of the permittivity, as shown in 3(c). It is important to note that in the dispersion and trajectories plots we only show the lower and upper polariton trajectories, no other modes are shown for the sake of clarity, but they also exist. with high f, and for hBN. The uncoupled modes are shown in grey dashed lines and the location of the excitonic resonance is shown in a yellow star. These hybridized modes are much harder to follow, in comparison with the waveguide modes because the uncoupled optical mode has a non-trivial trajectory (showed in dashed lines) with high losses (Im(ω). Thus, when it is coupled to the exciton it is attracted by it with a non-trivial trajectory. An additional issue with this modes is the cut-off induced by the light-line, such that after a certain k x the modes become guided. Fig. S5 depicts the analogue information as in the Fig. 3 in the main text, but for TM modes. In this case, the dispersion and the trajectories of the poles are given by the self-hybridization of the J-aggregate slab with the TM modes both below and above the light line. The 1 st FP mode is so weakly coupled that the trajectories barely change from the uncoupled ones (Fig. S5(e)), as expected there is no splitting noticeable in the dispersion (Fig. S5(a)). Nevertheless, TM 1 is strongly coupled, the coupling strength is part of the summary shown in Fig. 6 in the main text. On the other hand, when f ω 2 P = 3.116 eV 2 , there is splittig in the FP mode, thus the optical mode is clearly pulled towards the exciton in the trajectories plot (S5(f)). The splitting with TM 1 is so big that the mode at zero detuning for the upper polariton is cut by the light-line.

III. 1D: LONG CYLINDER
This section contains the calculations corresponding to infinite perovskite CsPbCl 3 cylinders with r = 250 nm, but showing the self-hybridization with the TEM and TE modes, in analogy to in Figure 4 in the main text. The figures do not show all the modes in the plotted energy spectra. The only radiating mode for TEM in Fig. S6 the plot is the 4 th mode, which shows the Rabi splitting. TEM has many guided modes that couple to the resonance, the only modes shown in the plot are the HE11 mode (the first guided one) and the 5th and 6th modes, which are the last modes crossing the resonance, ω 0 = 3 eV due to the light-line cut. As we can see in Fig. S6(c), all the modes below the light-line follow the same trajectories but at different values of k x . The Rabi splitting above the light-line is shown in an orange arrow and below it with a purple arrow. Both values of the Rabi splitting were considered for the summary plot in Fig. 6(b) in the main text. Fig. S7 shows the information (same material and dimensions) but for TE modes. In this case the radiative modes closest to the resonance are the 2 nd (•) and 3 rd (♦) modes. Both modes are highly detuned, thus we cannot see the Rabi splitting above the light-line for this radius. The detuning is also visible in the trajectories (Fig. S7(b)), where the modes are barely modified with respect to the uncoupled ones (dashed lines and star). On the other hand, below the light-line there is a clear Rabi splitting marked with a purple arrow for both TE 01 , , and TE 02 , ▽.

IV. 0D: POLARITONS IN SPHERES
This section shows the dispersion and the trajectories of the TM 1,1 mode of water spheres. Since the loss of the optical mode is too large (in comparison to the resonance loss γ), the optical mode and the resonance are only weakly coupled. This causes an crossing in the dispersion Re(ω) vs radius (a similar behavior was observed in Fig. 2(a)). On the other hand, in the complex-ω plane we observe just a small deviation of the trajectory of the optical mode and a very small circular trajectory around the ω 0 − iγ/2 (similar to the trajectory in Fig.2(d)). The dispersion in the case of spheres is given by different sizes of spheres, thus the colormap in this case shows the radius of the sphere in µm. V. LIMIT OF RABI SPLITTING Table I contains the calculated values of the bulk Rabi splitting, 2g 0 = ω P f /ε ∞ , for the studied permittivities. Such values were used to normalize the obtained Rabi splitting of the differ-ent structures of Fig. 6 in the main text. Due to the difference in materials, the structures (spheres, cylinders and slabs) have different dimensions.
For the 2D case, the Rabi splitting for perovskites was obtained for slabs with a thickness of L = 106 nm; for water we used L = 1.1 µm; and as mentioned in the text before, for J-aggregates we used L = 200 nm and for hBN L = 1.75 µm.
For the 1D case, the Rabi splitting for water was obtained from cylinders with a radius of r = 2 µm; for J-aggregates we used r = 400 nm; and for pervoskites we have mentioned already that r = 250 nm.
For the 0D case, spheres of different radius were used for obtaining the Rabi splitting of different modes. The smallest spheres are coupled to TE 1,1 . Thus this is the minimal radius for spheres sustain polaritons. In the case of perovskites that radius was r = 98 nm; for J-agregates it was r = 180 nm and for water we have mentioned before that r = 1 µm.