Sub-picosecond thermalization dynamics in condensation of strongly coupled lattice plasmons

Bosonic condensates offer exciting prospects for studies of non-equilibrium quantum dynamics. Understanding the dynamics is particularly challenging in the sub-picosecond timescales typical for room temperature luminous driven-dissipative condensates. Here we combine a lattice of plasmonic nanoparticles with dye molecule solution at the strong coupling regime, and pump the molecules optically. The emitted light reveals three distinct regimes: one-dimensional lasing, incomplete stimulated thermalization, and two-dimensional multimode condensation. The condensate is achieved by matching the thermalization rate with the lattice size and occurs only for pump pulse durations below a critical value. Our results give access to control and monitoring of thermalization processes and condensate formation at sub-picosecond timescale.


Introduction
Ideal Bose-Einstein condensation means accumulation of macroscopic population to a single ground state in an equilibrium system, with emergence of long-range order. Bosonic condensation in non-or quasi-equilibrium and driven-dissipative systems extends this concept and offers plentiful new phenomena, such as loss of algebraically decaying phase order (1,2), generalized Bose-Einstein condensation (BEC) into multiple states (3), rich phase diagrams of lasing, condensation and superradiance phenomena (4,5), and quantum simulation of the XY model that is at the heart of many optimization problems (6). Such condensates may also be a powerful system to explore dynamical quantum phase transitions (7). Each presently available condensate system offers different advantages and limitations concerning studies of non-and quasiequilibrium dynamics. The ability to tune interactions precisely over a wide range is the major advantage of ultracold gases (8,9). Polariton condensates (10)(11)(12)(13)(14)(15)(16)(17)(18) offer high critical temperatures compared to ultracold gases. In photon condensates (19,20), the thermal bath is easily controlled (21,22). Recently periodic two-dimensional (2D) arrays of metal nanoparticles, socalled plasmonic lattices or crystals (23), have emerged as a multifaceted platform for room temperature lasing and condensation at weak (24)(25)(26)(27)(28) and strong coupling (29,30) regimes.
In plasmonic lattices, the lattice geometry and periodicity, the size and shape of the nanoparticles, and the overall size of the lattice can be controlled with nanometer accuracy and independent of one another. The energy where condensation or lasing occurs is given by the band edge energy that depends on the period of the array. Remarkably the band edge energy and the dispersion are extremely constant over large lattices (accuracy 0.1% (28)). In semiconductor polariton condensates, disorder in the samples often leads to traps and fragmentation (15), or condensates may be trapped by geometry (16,19). Thus plasmonic lattices offer a feature complementary to other condensate systems, namely that propagation of excitations over the lattice can be used for monitoring time-evolution of such processes as thermalization: each position in the array can be related to time via the group velocity, and there are no spurious effects due to non-uniformity of the sample. Spatially resolved luminescence was utilized in this way in the first observation of a BEC in a plasmonic lattice (28).
Here we show that formation of a condensate with a pronounced thermal distribution is pos-sible at a 200 fs timescale and attribute this strikingly fast thermalization to partially coherent dynamics due to stimulated processes and strong coupling. We observe a unique double threshold phenomenon where one-dimensional (1D) lasing occurs for lower pump fluences and twodimensional (2D) multimode condensation, associated with thermalization, at higher fluences.
The transition between lasing and condensation shown in our work is different from previous condensates (16,28,(31)(32)(33): it relies on matching the system size, propagation of excitations, and the thermalization dynamics. Importantly we find a peculiar intermediate regime showing thermalization features but no macroscopic population at the lowest energy states. This regime allows us to reveal the stimulated nature of the thermalization process by the behavior of the luminescence in lattices of different sizes. As a direct evidence of the ultrafast character of the thermalization and condensation process, we show that it occurs only for pump pulse durations below a critical value of 100−250 fs.

System
Our system consists of cylindrical gold nanoparticles in a rectangular lattice overlaid with a solution of organic dye molecule IR-792 (see details in Section S1). The lattice supports dispersive modes, so-called surface lattice resonances (SLRs), which are hybrid modes composed of localized surface plasmon resonances at the nanoparticles and the diffracted orders of the periodic structure (23,34). The electric field of the SLR modes is confined to the lattice plane in which the SLR excitations can propagate. An SLR excitation can be considered a bosonic quasiparticle that consists (mostly) of a photon and of collective electron oscillation in individual metal particles. The dye molecules are pumped optically, and they can emit to the SLR modes.
We excite the sample with laser pulses at 1 kHz repetition rate and central wavelength of 800 nm, and resolve the luminescence spectrally as a function of angle and spatial position on the array, Fig. 1a (for details see Section S2). The SLR modes are classified to transverse magnetic (TM) or transverse electric (TE) depending on the polarization and propagation direction, as defined in Fig. 1b. The measured dispersions are displayed in Fig. 1c-e. In the presence of dye molecules, the SLR dispersion shifts downwards in energy with respect to the initial diffracted order crossing. Moreover the TE modes begin to bend when approaching the molecular absorption line at 1.53 eV. These observations indicate strong coupling between the SLR and molecular excitations (35); a coupled modes fit to the data gives a Rabi splitting of 164 meV (larger than the average line width of the molecule (150 meV) and SLR (10 meV)) and an exciton part of 23% at k = 0. In the following, we refer to the hybrids of the SLR mode excitations and molecular excitons as polaritons, for brevity. To our best knowledge, this is the first time that high molecule concentration in a liquid gain solution is implemented in plasmonic lattice systems, cf. (26,29). The plasmonic lattice, optimized for creating the condensate, has a particle diameter of 100 nm and height of 50 nm, the period in yand x-direction of p y = 570 nm and p x = 620 nm, dye concentration of 80 mM, and a lattice size of 100×100 µm 2 . The period p y is varied between 520 and 590 nm, and the lattice size between 40×40 and 200×200 µm 2 .

Results
We study the luminescence properties of the plasmonic lattice as a function of pump fluence, that is, the energy per unit area per excitation pulse, and find a prominent double threshold behaviour. The system is excited with an x-polarized 50 fs laser pulse that has a flat intensity profile and a size larger than the lattice. Excited polariton modes continuously leak through radiative loss, and therefore the observed luminescence intensity is directly proportional to the population of the polaritons. We record real space and momentum (k-)space intensity distributions and the corresponding spectra, the photon energy is E = hc/λ 0 and the in-plane wave vector k x,y = 2π/λ 0 sin(θ x,y ).
The sample luminescence as a function of pump fluence is presented in Fig. 2. The total luminescence intensity reveals two non-linear thresholds and a linear intermediate regime, shown in Fig. 2a. The line spectra, shown as insets, are obtained by integrating the real space spectra along the y-axis and unveil the population of polaritons as a function of energy. At the first threshold, Fig. 2b-c, lasing typical for nanoparticle arrays (24,27,29) is observed throughout the array. Increasing the pump fluence beyond the first threshold, Fig. 2d-e, the luminescence becomes more intense in the central part compared to the top and bottom parts of the array.
Moreover luminescence at the center takes place at a lower energy than at regions closer to the edges. We interpret this red shift as a signature of thermalizing population of polaritons that propagate along the array in +y and −y directions, discussed below. At the second threshold, Fig. 2f-g, the system undergoes a transition into a condensate: the real space intensity distribution shows uniform luminescence in the central part of the array, and the line spectrum (Fig. 2a, top inset) has a narrow peak at the band edge and a long thermalized tail at higher energies.
A fit of the tail to the Maxwell-Boltzmann distribution (dashed line) gives the room temperature, 333±12 K (see Section S3 for more information). The spectrometer counts per emitted condensate pulse correspond to a photon number of ∼10 9 (see Section S4), which is roughly 10 5 times more than in the first BEC in a plasmonic lattice (28). This tremendous improvement has increased the signal-to-noise ratio so that we can now observe a prominent thermal tail (it is likely to be even longer but the data is cut due to filtering out the pump pulse). Such an upgrade of luminescence intensity is crucial for future fundamental studies and applications of this type of condensate. To produce a condensate, the periodicity must be tuned with respect to the thermalization rate and the array size (Section S5).
Three distinct regimes are also observed in the k-space intensity distributions. Fig. 3 presents the k-space images and spectra for the same sample and pump parameters as in Fig. 2. Fig. 3ac shows that lasing spreads in the TM mode to large k (i.e., large emission angles), whereas thermalization of the polaritons occurs mostly along the TE mode, see Fig. 3d-f. At the condensation threshold, Fig. 3g-i, we observe confinement in both k x and k y , implying that the condensate has a 2D nature, in contrast to the lasing regime where confinement is observed only in k y . Fig. 3j shows line spectra obtained by integrating along k y from TE mode crosscuts ( Fig. 3c,f,i); the spectrometer slit width of 500 µm corresponds to ±1.3 • around θ x = 0 in the 2D k-space images. Intriguingly, at the condensation threshold, we see multiple modes highly occupied at k y = 0. The line spectrum at 3.49 mJ/cm 2 shows three narrow peaks followed by a thermalized tail. The full-width at half-maximum (FWHM) of the highest peak is 3.3 meV, significantly narrower than the bare SLR mode (10 meV). The spacing of the multiple peaks decreases towards lower energy, ruling out the possibility of a trivial Fabry-Pérot interference. We have also observed that the spacing does not depend on the lattice size or periodicity. Based on T-matrix simulations, we attribute these modes to the finite size of the nanoparticle (Section S6).
This highlights the role of the nanoparticles beyond providing a periodic structure.
Lasing and condensation transitions are expected to result in increased spatial and temporal coherence of the emitted light. Increase in temporal coherence was evidenced as the narrowing of spectral line width (Fig. 3j). To study the spatial coherence, we have performed a Michelson interferometer experiment as a function of pump fluence. In the Michelson interferometer, the real space image is split into two, one image is inverted and combined with the other one at the camera pixel array. The contrast of the observed interference fringes is extracted with a Fourier analysis of the spatial frequencies in the interfered image to exclude any artifacts produced by intensity variations in a single non-interfered image (Section S7).
In Fig. 4, the fringe contrast (proportional to the first-order correlation function g (1) ) is shown as a function of pump fluence in both yand x-directions of the lattice. High spatial coherence occurs in the y-direction over the array in both the lasing and condensation regimes ( Fig. 4b-c). At the intermediate regime, the spatial coherence decreases, in agreement with the observation that the thermalizing population dominates the luminescence signal. In the x-direction, spatial coherence is lower than in the y-direction over the whole pump range.
However, the condensate clearly exhibits high spatial coherence also in x-direction, in contrast to lasing, which shows separated emission stripes in the real space image (see Fig. 4d and Fig. S5). The spatial coherence measurements are in line with the observations from the 2D k-space images (Fig. 3), where 1) lasing exhibits confined luminescence along k y (the direction of feedback) but spreads along k x , 2) condensation shows 2D k-space confinement.
The lasing and condensation take place at energies 1.397 and 1.403 eV, respectively, lower than the band edge energy of the system without molecules, 1.423 eV. However the energies are blue-shifted from the lower polariton energy 1.373 eV (band edge in reflection) or 1.382 eV (fitted lower polariton branch; Fig. 1d), which are experimentally obtained by a reflection measurement and coupled modes model (35). Since the whole dispersion gradually blue shifts as a function of pump fluence, it may be associated to degradation of strong coupling instead of polariton-polariton interactions. Based on the band-edge locations (1.423 -1.403) eV / (1.423 -1.373) eV, the coupling in the condensation regime has decreased to ∼40% of the case without pumping. Note that the observed double-threshold behaviour is different from semiconductor microcavity polariton condensates where condensation has a lower threshold than photonic lasing associated with loss of strong coupling (16,31). However, careful comparison of the distance between the array edge and the location where the red shift begins reveals that for all arrays, for the given pump fluence, the distance is the same (∼ 25 µm; see Fig. 5a,c,e for definition of the distance).
We explain this distance by stimulated emission pulse build-up time: the time between the maxima of the population inversion and the output pulse as defined in the rate equation simulation in Fig. 5g (see Ref. (36) and Section S8 for description of the model). Pulse buildup is a well-known phenomenon in Q-switched lasers (37). In our system, the pump pulse excites the molecules, and the polaritons begin to propagate when the first spontaneous photons populate the modes. The modes then gather gain while propagating, and therefore the peak of the stimulated emission pulse appears after a certain distance travelled along the array. This distance is seen as the dark zones in real space measurements, and it corresponds to the pulse build-up time. Note that this pulse is different to a lasing or condensate pulse since the line width of the luminescence is large and spatial coherence is small. By summing up such spatial intensity profiles of the thermalizing pulses at every point along the lattice (Fig. 5h), we can reconstruct the real space intensity distributions (insets of Fig. 5a,c,e). The dark zones appear because the edges do not receive propagating excitations from outside the lattice; the dark zones at the edges have approximately half the intensity of the central part. In the small lattice, intensity at the edges is similar to that of the larger lattices but in the center it is only half of that. Moreover the wavy interference patterns in the central part (Fig. 5c,e) indeed only appear for arrays larger than 40×40 µm 2 , where there are counter-propagating pulses that can interfere.
We found that the width of the dark zones depends on pump fluence as predicted by the rateequation simulation and the theory of Q-switched lasers, namely the build-up time is inversely proportional to the pump fluence. Fig. 5i shows that the dark-zone width follows the inverse of the pump fluence until it saturates at around 3 mJ/cm 2 to a value below 20 µm (∼100 fs). The inset in Fig. 5i shows the pulse build-up time extracted from our rate-equation simulation, and it displays a similar ∼1/P dependence.
We attribute the red shift to a thermalization process. At the intermediate regime of pump fluences, however, a thermal distribution is not reached before the population decays. For higher fluences, a condensate peak and a tail with Maxwell-Boltzmann distributed population emerges.
Emission-absorption cycles with dye molecules are known to provide thermalization of photons (19,32,38) and lattice plasmons (28), at the weak coupling regime, to the temperature of the vibrational degrees of freedom if several (39) (or at least one (20)) cycles take place within the cavity lifetime. In some studies of organic polariton condensates, condensation is attributed to loss of a vibrational quantum (14,29,40,41). The SLR mode lifetimes are only about hundred femtoseconds, so it is unprecedented that we observe a prominent thermal distribution with dynamic range over one decade. We propose that this is possible because the thermalization is highly coherent due to the observed stimulated processes and due to strong light-matter coupling. For a single molecule with one vibrational state coupled to a few light modes, one can theoretically describe how a process that is coherent except for vibrational losses leads to rapid red shift of emission (Section S9). Describing quantum dynamics of molecules with thermal vibrational state populations, interacting strongly with multiple modes of light at the multi-photon regime, is beyond the current state-of-the-art theory (42)(43)(44)(45). Such an advanced description will be needed to rigorously describe our observations.

Effect of pulse duration
The spatial measurements have enabled an astounding, yet indirect, way to look into the dy- Interestingly this is similar to or smaller than the time (∼250 fs) in which the polaritons propagate from the edges to the center in the 100×100 µm 2 array. For a 250 fs pulse duration, some narrow peaks occur at the band edge but the thermalization in the time-integrated signal remains incomplete (too much population in the higher energy states), which can be because of molecules are brought to the excited state also during the thermalization process. With the longest excitation pulses (>350 fs), it seems the instantaneous population inversion does not reach a high enough value to even start the thermalization process as it competes from the same gain with the lasing, already triggered at the first threshold (see also Fig. 6). Sensitivity to the excitation pulse duration highlights the ultrafast nature of plasmonic systems and endorses the sub-picosecond dynamics of the thermalization process. Besides the critical pulse duration, the condensation requires that the thermalization time matches the propagation time so that the polaritons have red shifted to the band edge energy while still having large population density, which can be achieved with an optimal balance between the dye concentration, lattice size, and periodicity.

Conclusions
Fundamental questions on the dynamics of Bose-Einstein condensation in driven-dissipative systems are still largely open, despite years of research (1,2,15,16). What is the nature of the energy relaxation and thermalization processes, how does the condensate form, and what are its quantum statistical and long-range coherence properties? These questions become challenging to answer for room temperature condensates as higher energy scales imply faster dynamics.
Here we have shown that plasmonic lattices offer an impressive level of access and control to the sub-picosecond dynamics of condensate formation via propagation of excitations and the finite system size.
We experimentally demonstrated that a bosonic condensate with a clear-cut thermal excited state population can form in a timescale of a few hundreds of femtoseconds. We propose that this extraordinary speed of thermalization is possible because the process is partially coherent due to strong light-matter coupling and stimulated emission. Strong light-matter coupling at the weak excitation limit was indicated by our reflection measurements. By varying the lattice size, we revealed the stimulated nature of the thermalization process. While strong light-matter coupling at the multi-photon regime is described by the Dicke model (46), to explain both the red shift and the thermal distribution we observed, one would need to involve vibrational degrees of freedom that are (strongly) coupled to the electronic ones as well as to a thermal bath. Work toward surmounting such theoretical challenges has already begun since systems where light and molecular electronic and vibrational states are strongly coupled have promise for lasing and condensation, energy transfer, and even modification of chemical reactions (47,48). We have shown here that plasmonic lattices offer a powerful platform for studies of such ultrafast light-matter interaction phenomena. The shape, size, and material of the nanoparticles can be accurately tuned, as well as the lattice geometry, composition of the unit cell, and the lattice size.
This provides a large, controlled parameter space vital for testing and (dis)qualifying theoretical predictions. In particular, the dynamics can be accessed in complementary ways: through conventional time-domain techniques as well as indirectly via the propagation of excitations in the lattice.

Acknowledgments
We thank Janne Askola for his help in intensity calibration of the spectrometer. We thank 0.83 mJ/cm 2 1.47 mJ/cm 2 3.49 mJ/cm 2           The gold nanoparticle arrays are fabricated with electron beam lithography on glass substrates where 1 nm of titanium is used as an adhesion layer (see SEM image in Fig. S1). The nominal dimensions of the plasmonic lattice, optimized for the condensation experiment, are the following: a nanoparticle diameter of 100 nm and height of 50 nm, the period in yand x-direction of p y = 570 and p x = 620 nm, respectively, and a lattice size of 100×100 µm 2 . In reference measurements, the period p y is varied between 520 and 590 nm and the lattice size between 40×40 and 200×200 µm 2 .
Asymmetric periodicity separates the diffracted orders in the energy spectrum for the two orthogonal polarizations (e x and e y ), and the SLR dispersions are correspondingly separated, which simplifies the measured spectra. For x-polarized nanoparticles (as in our experiments), the TE and TM modes correspond to combinations of (e x , k y ) and (e x , k x ), respectively. Under pumping, which SLR mode is mainly excited is determined by the pump polarization as the molecules are excited more efficiently via the single particle resonance at the plasmonic hotspots of each nanoparticle (49). In the experiment with different periods, p x was kept always 50 nm larger than p y . In the experiment where lattice size was varied, however, the lattice period in x and y directions was the same (p x = p y = 570 nm). We found that asymmetric periodicity does not play a crucial role in forming the condensate but just simplifies the data analysis of the experimental results.
The dye molecule solution is index-matched to the glass substrate (n = 1.52), the solution being a mixture of 1:2 DMSO:Benzyl Alcohol. The solution is sealed inside a Press-to-Seal silicone isolator chamber (Sigma-Aldrich) between the glass substrate and superstrate. The dye solution has a thickness of ∼1 mm, which is very large compared to the extent of the SLR electric fields. Active region of the dye lies withing a few hundred nanometers from the lattice plane, shown experimentally in (27,36). Given by the excess of the dye molecules and the natural circulation of the fluid, there are always fresh dye molecules available for each excitation pulse, which makes the sample extremely robust and long-lasting. IR-792 perchlorate (C 42 H 49 ClN 2 O 4 S) was chosen as the dye molecule because it dissolves to the used solvent in very high concentrations, in contrast to many other dye molecules, e.g., IR-140 that has also been used by us (27,36) and others (24,26,49) as a gain medium in plasmonic nanoparticle array lasers.

Section S2. Transmission, reflection, and luminescence measurement setup
A schematic of our experimental setup is depicted in Fig. S2. The same setup is used for transmission, reflection and luminescence measurements with minor modifications. The spectrometer resolves the wavelength spectrum of light that goes through the entrance slit, and each pixel column in the 2D CCD camera corresponds to a free space wavelength, λ 0 , and each pixel row to a y-position at the slit. The y-position further corresponds to either an angle (k-space) or the y-position at the sample (real space). The photon energy is E = hc/λ 0 and in the case of angleresolved spectra (dispersions) the in-plane wave vector k x,y = k 0 sin(θ x,y ) = 2π/λ 0 sin(θ x,y ), where h is the Planck constant and c the speed of light in free space. Next, the three different experiment types are explained, starting with the luminescence measurement where the sample is optically excited with an external pump laser. The excitation pulse (or pump pulse) is generated by Coherent Astrella ultrafast Ti:Sapphire amplifier. The pulse has a central wavelength of 800 nm, and at the laser output, a duration of < 35 fs with a bandwidth of 30 nm. The pump pulse is guided thorough a beam splitter and mirrors, and finally to the mirror M 1 (see Fig. S2), which directs the pump pulse to the excitation path of the setup. We have a band-pass filter in the excitation path that is used in combination with a long-pass filter in the detection path to filter out the pump pulse in the measured luminescence spectra. The pump pulse is linearly polarized, and to filter only the horizontal component we have placed a linear polarizer after the band-pass filter. The pump fluence is controlled with a metal-coated continuously variable neutral-density filter wheel (ND wheel). The pump pulse is spatially cropped with an adjustable iris and the iris is imaged onto the sample with a help of lens L1 and the microscope objective. The inverted design enables exciting the sample at normal incidence, which is crucial for simultaneous excitation of the dye molecules over the sample. Excitation at normal incidence also prevents any asymmetry in the spatial excitation of the molecules around the nanoparticles with respect to lattice plane.
The inverted pumping scheme and accurate optical alignment were essential for repeatable and precise condensate formation.
In the detection path, we have the long-pass filter and optionally a linear polarizer. An iris or pinhole acts as a spatial filter at the 1st image plane to restrict the imaged area at the sample.
The 1st image plane is relayed onto the real-space camera (1st Cam.). In the k-space measurements, the back-focal plane of the objective (Fourier plane; containing the angular information of the collected light) is imaged to the 2nd Cam., that acts as a 2D k-space camera, and to the spectrometer slit, with the tube lens and a k-space lens. For the real space measurements, the beam-splitter before the k-space lens is replaced with an additional real-space lens to produce the 2nd image plane of the sample to 2nd Cam. and onto the spectrometer slit. The spectrometer slit selects a vertical slice either of the 2D k-space image or the real space image. In the luminescence measurements, we use a slit width of 500 µm. In the k-space, it corresponds to  from a halogen lamp. The lattice modes are visible as transmission minima (extinction maxima) in the angle-resolved spectrum. When a thick layer of high-concentration dye molecule solution is added in the chamber on top of the nanoparticle array, the transmission measurement is not applicable due to a complete absorption by the molecules. To access the dispersion of the lattice modes in this case, we use the setup in reflection mode by utilizing the same inverted design as in the luminescence measurement. The halogen lamp is inserted before the iris in the excitation path, that is imaged onto the sample, and the dispersion of the lattice modes is revealed by reflection (scattering) maxima in the collected angle-resolved spectrum.
The luminescence measurement as a function of pump fluence is automated with LabView.
Predefined fluence steps are measured such that for each step: 1) the calibrated ND filter wheel is set to a correct position, 2) the shutter is opened, 3) the image is acquired with spectrometer and optionally with 1st and 2nd Cam., and 4) the shutter is closed. Integration time of the spectrometer is automatically adjusted during the measurement to avoid saturation at highly non-linear threshold regimes. The pump pulse duration is measured with a commercial autocor-relator (APE pulseCheck 50). In the pulse duration measurement, the pump pulse goes through the same optics as in the actual experiments. The pulse duration is changed by adjusting the stretcher-compressor of the external pump laser.

Section S3. Fits to the Maxwell-Boltzmann distribution
We fit the thermalized tail in the measured population distributions to Maxwell-Boltzmann distribution (manuscript Fig. 2 and Fig. S8). The fit function is given by 3.9 mJ/cm 2 , the fits to the the Maxwell-Boltzmann distribution give temperatures of 291±7 K and 262±9 K, respectively. For pump fluences above ∼4 mJ/cm 2 , the linear slope is distorted and the condensate degrades. This is evident also as a decrease of the spatial coherence (manuscript Fig. 4a) and an increase of the FWHM of the spectral maximum (Fig. S8c).

Section S4. Estimation of photon number in the condensate
The photon number in the condensate is estimated from the measured luminescence intensity.
A strongly attenuated beam from the external pump laser (800 nm, 1 kHz) is directed to the spectrometer slit, and the total counts given by the spectrometer CCD camera (Princeton Instruments PIXIS 400F) is compared to the average power measured with a power meter (Ophir Vega). The measured average power of 167 nW corresponds to 6.7 · 10 8 photons/pulse whereas the total counts in the CDD camera are 8.4 · 10 6 , leading to a conversion factor of ∼80 photons/count. In the condensation regime, the total counts per pulse is about 3 · 10 6 (manuscript Fig. 2a). The collection optics including the beam splitters reduce the signal roughly by a factor of 2.5, and as the slit width of 500 µm corresponds to 27 µm at the sample, we collect luminescence from an area that is about 1/4 of the 100 µm wide nanoparticle array. Finally, the sample is assumed to radiate equally to both sides, so the actual photon number per emitted condensate pulse becomes: n ph ≈ 2.5 × 4 × 2 × 80 × 3 · 10 6 = 4.8 · 10 9 .

Section S5. Band-edge energy
We demonstrate the effect of the band edge location in Fig. S3 with two examples where the band edge is tuned either below or above the optimal energy of 1.40 eV. When the periodicity is set small so that the band edge of the x-polarized TE mode is at high energy ( Fig. S3a-c), the propagation and red shift occurs along the lower dispersion branch of the TE mode, i.e., the red shift is not halted at the band edge. In contrast, when periodicity is set large so that the band edge resides at low energy ( Fig. S3d-f), we observe the red shift of polaritons propagating along the upper dispersion branch of the TE mode but no condensation. The redshifting population simply does not reach the band edge, which is at too low energy. The stripes in the real space spectra, Fig. S3b,e (also visible in manuscript Fig 2d,e and Fig. 4c-f), arise from standing waves caused by counter-propagating modes. We found that the wavelength of the intensity oscillations is λ RS (E) = π/k(E). Comparing the k-space and real space spectra shows in Fig. S3a-b that the oscillations are denser at lower energies because the lower energy corresponds to larger k. In contrast in Fig. S3d-e, the oscillations are sparser at lower energies because the lower energy corresponds to smaller k.
Interestingly while in previous studies using a dye molecule bath for photons (19,32) or lattice plasmons (28) condensation required matching the lowest state energy with the energy where the molecule absorption vanishes, here that condition is not needed but the condensate formation is controlled by the interplay of the lattice size, periodicity and thermalization speed.
Counts per pulse (#) Counts per pulse (#) Figure S3: Effect of the band-edge energy. Left column: k-space spectra in the TE mode direction.
Middle column: Real space spectra. Right column: Line spectra obtained by integrating over the real space spectra in y-direction (between white lines). Results are shown for two lattice periods: (a-c)

Section S6. T -matrix simulation
To unveil the origin of the multiple modes, we have performed multiple-scattering T -matrix simulations of an infinite array of spherical nanoparticles, with the periods in x and y directions corresponding to our system.
In the T -matrix approach, the scattering properties of a single nanoparticle are first described in terms of vector spherical wave functions (VSWFs), giving the T -matrix of the particle at a given frequency. For a spherical particle, the nontrivial elements of the T -matrix are given by the well-known Lorentz-Mie solution (50). Next, the interactions between nanoparticles at different positions are expressed in terms of translation operators between the VSWFs with different origins (51). Applying the Bloch boundary conditions, the electromagnetic response of a periodic nanoparticle array can be described with a matrix equation of the form where T = T (ω) is the single particle T -matrix, W (ω, k) is a lattice sum of the VSWF translation operators, a is a vector containing the coefficients of VSWFs scattered from a nanoparticle, and p ext is a vector of incoming VSWF coefficients, describing the external fields driving the array. Lattice modes are defined as the solutions of Eq. (2) with the right hand side set to zero (physically this means the waves propagate without the need of external driving) and exist only for such (ω, k) pairs for which the matrix (I − T W ) is singular (giving the dispersion relation of the array), which is equivalent to some singular value of the matrix being equal to zero. In practice, the particles are lossy, hence the frequency ω needs to be complex if the wave vector k is real. We also exploit the symmetries of the system and evaluate the equation (2)   radii. For the smallest nanoparticle, the singular value minima (which correspond to the lattice modes) of all irreducible representations are concentrated at one energy. As the particle size is increased, some of the minima (hence also the modes) start to separate from each other in energy. The results show that multiple degenerate modes exist at k = 0 for small particles, but larger particle size lifts this degeneracy producing closely spaced modes with non-uniform spacing, see Fig. S4. We conclude that the multi-peak structure of the condensate can thus originate from the complexity of finite-sized nanoparticles.

Section S7. Spatial coherence measurements
Spatial coherence of the sample luminescence is measured with a Michelson interferometer.
The real space image is split into two arms and the image in one of the arms is inverted in vertical direction with a hollow roof retro-reflector. Then the two images are combined again with a beam splitter and overlapped at the camera pixel array, simultaneously. With this design the spatial coherence (first-order correlation g (1) ) can be measured separately along both xand y-axis of the plasmonic lattice. The retro-reflector always inverts the image vertically, with respect to y = 0 in laboratory reference frame, but the sample and the pump polarization can be rotated 90 • to measure g (1) (−x, x) instead of g (1) (−y, y) in the lattice coordinates. The first-order correlation function describing the degree of spatial coherence is given by where E(y) is the electric field at point y. The first-order correlation function relates to the interference fringe contrast C as follows: C(y, −y) = 2 I(y)I(−y) I(y) + I(−y) g (1) (y, −y), where I(y) is the luminescence intensity at point y of the lattice. The fringe contrast in the interfered images is extracted with a Fourier analysis as explained below.
First we need to determine the period of the interference fringes arising due to coherence of E(y) and E(−y), and set a Fourier filter for spatial frequencies accordingly. After that, the image data is gone through column by column, inside the region of interest, and a Matlab inbuilt Fast-Fourier Transform algorithm is performed to each pixel column at a time. In the spatial frequency spectrum, we find the peak value inside a predefined frequency bandwidth and compare that to the noise floor. The peak value needs to be above a chosen threshold value. If the peak is above the threshold, the rms-sum of the frequency components within the predefined bandwidth is compared to the background level (DC value of the frequency spectrum). Finally the contrast at certain pump fluence is taken as the mean value of the contrast along the columns inside the region of interest. The threshold value for a peak-acceptance level is adjusted so that the mean contrast value found in the reference cases (incoherently summed real space image) stays below 5%. Fig. S6 shows an explanatory example when the method is not applied for each pixel column separately but the intensity along y-axis averaged over x (from Fig. S5b), and the intensity along x-axis averaged over y (from Fig. S5d).  For that, Fourier analysis of spatial frequencies is necessary to extracting the fringe contrast caused by spatial coherence. Here P th means the threshold pump fluence for condensation and E refers to the pump polarization. The real space images are recorded for single pump pulses.  Fig. S5b), and the right panel (b,d,f) shows the intensity along x-axis of the lattice, averaged over y (Fig. S5d). In (c-f), the blue curves show the original signal and the red curves the filtered signal.

Section S8. Rate-equation simulation of a stimulated-emission pulse
We use a standard four-level model for the gain medium to simulate the stimulated emission pulse when the four-level system is originally in its ground state, and excited with a 50 fs pump pulse. The levels are labeled as follows: the pump excites the system from the level 0 to 3, there is a non-radiative decay from 3 to 2, and emission to the cavity mode is from 2 to 1. The model shows the same temporal evolution after the pump pulse as a gain-switched laser, or a Qswitched laser after the Q-switch is opened (37). The transition lifetimes used for the four-level gain medium are: τ 32 = τ 10 = 50 fs and τ 21 = τ 20 = 500 ps, which are similar to those used in the literature for organic dye molecules (24,27,36). The spontaneous emission coupling factor (β-factor) is set to β = 0.001 and the cavity lifetime to τ cav = 100 fs (corresponding to a typical lifetime of an SLR mode). The model is defined with the following coupled rate-equations, as in Ref. (36): where the populations of each level are denoted with N i and the transition lifetimes with τ i .
Here, n ph is the photon number in the mode (in our case the number of polaritons). Parameter r is the pump rate proportional to the pump intensity that has a Gaussian temporal shape.
In the model, the threshold value for population inversion is defined by comparing the gain and loss terms for the photon number in Eq. (5). The optical gain must overcome the loss, and at the threshold they are equal Section S9. Description of the quantum model We have studied the thermalization mechanism qualitatively with a microscopic quantum model including multiple cavity modes coupled to a single two-level system that is coupled to a shifted harmonic oscillator that describes the rotational-vibrational degrees of freedom within a molecule. The results of the model are presented in Fig. S7.
The system is described by the Holstein-Tavis-Cummings model (43,(53)(54)(55)(56) with the Hamil- Here a † i is the bosonic creation operator of the cavity mode of index i, σ z and σ ± are the Pauli operators describing the two-level structure of the molecule and b † is the bosonic creation operator corresponding to the vibrational mode of the molecule. Furthermore, ω k is the energy of the cavity mode of index i, ω v is the energy of the vibrational mode, ω m is the energy of the two-level system, g is the coupling between cavity modes and the molecule, and S is the Huang-Rhys parameter. Rotating wave approximation has been used in the simulation, assuming that the coupling |g| is significantly smaller than the molecule and cavity mode frequencies. This is true with the parameters used in the simulation: ω i = 1. Dissipations of an open quantum system are taken into account in the Lindblad formalism, which yields the master equation (57): where Section S10. Different pump pulse durations We studied the condensation phenomenon as a function of pump pulse duration and found that thermalization and condensation happens only for sub-250 fs pump pulses. Comparison of 50 fs and 500 fs pump pulses is presented in Fig. S8. It is evident that the longer excitation pulse results in only one (lasing) threshold. The distributions at around the threshold (blue and red curves in Fig. S8a,e are similar with both pulse durations, but at the higher fluences (yellow, purple, and green curves) the distributions are very different. With a 500 fs pulse, the population does not reach a thermal Maxwell-Boltzmann distribution, and no narrow peaks appear at the band edge. Besides the luminescence intensity, the different threshold behaviour is clearly visible in the FWHM curves (Fig. S8c,g). The FWHM is significantly decreased with both pulse durations at the first (lasing) threshold but only in the 50-fs case, the FWHM is decreased even further at the second (condensation) threshold. Note that at intermediate pump fluences, the 50 fs pulse shows a sharp increase of the FWHM because the maximum of the line spectra is found at higher energies (see Fig. S8d) as the thermalizing population dominates the signal.
Examples of real space and k-space images and spectra for the 500 fs are shown in Fig. S9, for pump fluences corresponding to the (a,c) lasing threshold and (b,d) beyond the condensation threshold of the 50 fs pump pulse. Importantly the real space and k-space images and spectra look nearly identical for both pump fluences for the 500 fs pulse duration, confirming that there indeed is no second threshold where the condensation would take place.