Bose-Einstein Condensation in a Plasmonic Lattice

TL;DR

This paper reports the first observation of Bose-Einstein condensation of surface plasmon polaritons in a nanoparticle lattice, achieved through ultrafast thermalization at room temperature, with potential for on-chip quantum technologies.

Contribution

It demonstrates a novel BEC of surface plasmon polaritons in a lattice, with ultrafast thermalization and a tunable crossover to lasing, advancing plasmonic quantum coherence research.

Findings

Observation of BEC in plasmonic lattice modes

Ultrafast thermalization within picoseconds

Crossover from BEC to lasing by band structure engineering

Abstract

Bose-Einstein condensation is a remarkable manifestation of quantum statistics and macroscopic quantum coherence. Superconductivity and superfluidity have their origin in Bose-Einstein condensation. Ultracold quantum gases have provided condensates close to the original ideas of Bose and Einstein, while condensation of polaritons and magnons have introduced novel concepts of non-equilibrium condensation. Here, we demonstrate a Bose-Einstein condensate (BEC) of surface plasmon polaritons in lattice modes of a metal nanoparticle array. Interaction of the nanoscale-confined surface plasmons with a room-temperature bath of dye molecules enables thermalization and condensation in picoseconds. The ultrafast thermalization and condensation dynamics are revealed by an experiment that exploits thermalization under propagation and the open cavity character of the system. A crossover from BEC to usual lasing is realized by tailoring the band structure. This new condensate of surface plasmon lattice excitations has promise for future technologies due to its ultrafast, room-temperature and on-chip nature.

Equations

$$C(\mathbf{y},-\mathbf{y})$$

$$H_{\text{system}}=\sum _{\text{k}}\hslash {\omega }_{\text{k}}{a}_{\text{k}}^{†}{a}_{\text{k}}+\sum _{i=1}^N\left(\frac{\hslash {\omega }_m}{2}{\sigma }_i^z+\hslash \Omega b_i^{†}{b}_i-\hslash \Omega {\sigma }_i^z\sqrt{X}(b_i^{†}+b_i)\right)+\sum _{\text{k}}\sum _i^N\left(g_i{a}_{\text{k}}{\sigma }_i^{+}+g_i^{\ast }{a}_{\text{k}}^{†}{\sigma }_i^{-}\right).$$

$$\frac{d}{dt}{n}_{\text{k}}=-\kappa (\text{k})n_{\text{k}}+2N{g}^2\left[{\Gamma }^{+}(\text{k})(1+n_{\text{k}})p_e-{\Gamma }^{-}(\text{k})n_{\text{k}}(1-p_e)\right]$$

$$\frac{d}{dt}{p}_e=P(t)(1-p_e)-{\gamma }_m{p}_e+\sum _{\text{k}}2g^2\left[-{\Gamma }^{+}(\text{k})(1+n_{\text{k}})p_e+{\Gamma }^{-}(\text{k})n_{\text{k}}(1-p_e)\right].$$

$$I(\mathbf{r},{\mathbf{r}}^{\prime })=A/2⟨{\mathbf{E}}^{\ast }(\mathbf{r})\mathbf{E}({\mathbf{r}}^{\prime })⟩=A/2⟨|{\mathbf{E}}_1{|}^2+|{\mathbf{E}}_2{|}^2+2{\mathbf{E}}_1{\mathbf{E}}_2\cos \theta ⟩$$

$$g^{(1)}(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{⟨E^{\ast }(\mathbf{r})E({\mathbf{r}}^{\prime })⟩}{⟨E^{\ast }(\mathbf{r})⟩⟨⟨E({\mathbf{r}}^{\prime })⟩}.$$

$$C(\mathbf{y},-\mathbf{y})$$

$$E(\mathbf{k})$$

$$n_c$$

$${\lambda }_T(T)$$

Full text

Bosonic quantum statistics imply that below a certain critical temperature or above a critical density the occupation of excited states is strictly limited, and consequently, a macroscopic population of bosons accumulates on the ground state [1]. This phenomenon is known as Bose–Einstein condensation (BEC). Superconductivity of metals and high-temperature superconducting materials are understood as BEC of Cooper pairs [2, 3]. The BEC phenomenon is central in superfluidity of helium although the condensate constitutes a small fraction of the particles [4]. Textbook Bose–Einstein condensates with large condensate fractions and weak interactions were created with ultracold alkali atoms [5, 6, 7], and the fundamental connection between the superfluidity of Cooper pairs and the Bose–Einstein condensation was confirmed by experiments with ultracold Fermi gases [3]. While all these condensates allow essentially equilibrium description, as was the original one by Bose and Einstein, the phenomenology has expanded to non-equilibrium systems [8, 9, 10, 11, 12]. Hybrid particles of semiconductor excitons and cavity photons, called exciton-polaritons, have shown condensation and interaction effects [13, 14, 15, 16, 17, 18, 19], creating coherent light output that deviates from usual laser light. Magnons, that is, spin-wave excitations in magnetic materials [20, 21], and photons in microcavities [22, 23] form condensates as well. The most technologically groundbreaking manifestation of macroscopic population due to bosonic statistics has so far been laser light, which is a highly non-equilibrium state not thermalized to a temperature of any reservoir. As the BEC phenomenon has been observed only in a limited number of systems, new ones are needed for pushing the time, temperature and spatial scales where a BEC can exist, as well as for opening viable routes to technological applications of BEC.

Here we report the observation of BEC for bosonic quasiparticles that have not previously been shown to condense: excitations of surface lattice resonance (SLR) modes [24, 25, 26, 27, 28] in a metal nanoparticle array. SLRs originate from coupling between localized surface plasmon resonances of the nanoparticles and diffracted orders of an optical field incident on the periodic structure. The quasiparticle is mostly of photonic nature but also contains a matter part, namely electron plasma oscillations in the nanoparticles. This is different from the excitonic matter part in polariton condensates [9, 10, 12, 13, 14, 15, 16, 17, 18, 19], and distinguishes our system from the purely photonic condensates [22, 23]. The SLR excitations are essentially non-interacting, therefore we use weak coupling interaction with a molecular bath [22, 29, 23] to achieve thermalization. Due to the plasma oscillation component, losses fundamentally limit [30] the lifetimes of plasmonic systems [31, 32] typically to 10 fs – 1 ps. To achieve a condensate with Bose–Einstein statistics at room temperature, we have predicted [28, 33, 34] that thermalization timescales need to be 2–3 orders of magnitude shorter than in photon condensates [22] and similar to the fastest dynamics in polariton condensates [16, 18, 19]. Lasing has been observed in plasmonic lattices in the weak coupling [35, 36] and strong coupling (polariton) [37] regimes. To compare SLR lasing and BEC, we realize an experimental setting that provides direct access to the evolution of thermalization and condensate formation in the sub-picosecond scale.

The system

Our systems constitute of periodic arrays of gold nanoparticles on a glass substrate with inter-particle distances of 580–610 nm and array dimensions of 100 $\times$ 300 $(\mathrm{\mu m})^{2}$, see Figs. 1a-b, and Methods. Figure 1c shows the measured dispersion of the SLR mode. Due to coupling between the counter-propagating modes at ${\bf k}=0$ ($\Gamma$-point), a band gap opening typical for a grating is observed with a band edge in the upper dispersion branch. The curvature of the band bottom corresponds to an effective mass $m_{eff}\sim 10^{-37}$ kg ($10^{-7}$ electron masses), obtained by fitting to a parabola, Fig. 1c. The linewidth of the SLR mode near the band edge is 5.5 meV (3.8 nm) and the lifetime is 120 fs.

We use 50 mMol concentration molecule solution for which weak coupling between the SLR modes and the molecules is expected [38]. Figure 1e shows the absorption and emission profiles of the molecule (IR-792 perchlorate $\mathrm{C}_{42}\mathrm{H}_{49}\mathrm{ClN}_{2}\mathrm{O}_{4}\mathrm{S}$), and the point in energy ($\xi_{\mathrm{abs}}$), where absorption rate is effectively zero (i.e., small compared to the loss rate), is marked as a solid line. We excite the molecules with a 100 fs laser pulse and measure the far-field emission.

The SLR modes are confined via the near-field character of the nanoparticle plasmonic resonances as well as by interference effects due to the periodic structure. The majority of the field intensity lies close (within a few hundreds of nanometers) to the array surface. Thus the array is an open system that can be easily imaged and pumped from the top; changing the properties of the modes, such as quality factor, does not affect the direct access for pumping the molecules. We utilize this to monitor the spatial evolution of the excitations, see Figs. 1a-b.

Observation of Bose–Einstein Condensation (BEC)

We pump the dye molecules at the edge of the array and record the spectral evolution of SLR excitations propagating along the array, as sketched in Figs. 1a-b. We investigate three samples that are otherwise identical except that the inter-particle distance is varied to tune the energy of the SLR band edge a) above, b) below, and c) matching $\xi_{\mathrm{abs}}$. For intermediate cases, see Extended Data Fig. 1.

In Fig. 2a, we observe first emission around the emission maximum of the molecule and then a gradual decrease in energy (red-shift) when the population propagates along the array. When a photon is absorbed to one of the rovibrational levels of the molecule excited state, the excitation is likely to relax to some of the lower levels before a photon is emitted. Such recurrent emission and absorption cycles provide the red-shift. In the case of Fig. 2a, the band edge energy (dashed line) is well above $\xi_{\mathrm{abs}}$ (solid line). The molecules offer a bridge over the band gap, as photons are absorbed from the energies of the upper dispersion branch and emitted to the lower one. Eventually the red shift saturates to a broad spectral range at 1.34 eV, where the absorption has ceased (this is how we define $\xi_{\mathrm{abs}}\approx$ 1.34 eV).

The occupation of the dye molecule rovibrational states follow Maxwell–Boltzmann statistics at room temperature due the contact of the dye molecules with solvent molecules, thus the molecules provide a room-temperature thermal bath for the SLR excitations to exchange energy. Excitation by the pump may, however, drive the molecules out of equilibrium. Note that in our system, the SLR excitations propagate away from the initial excitation spot and there are always new molecules in thermal equilibrium available.

The spatial evolution over the array can be mapped to time via the group velocity of the SLR mode everywhere except a few meV range of energy near the band edge, where the group velocity is vanishingly small. Thus we can extract the thermalization time (or rate) from the slope of the propagating population. Absorption probabilities can be increased by number $N$ of the molecules, which leads to approximately linear dependence of the thermalization rate on $N$ [29]. We varied the molecular concentration and measured the thermalization rate, giving the expected linear dependence, Fig. 2e. In addition, the thermalization may be accelerated by the field hot spots near the nanoparticles (Fig. 1d), as Purcell enhancement factors of the order of hundreds have been reported for similar systems [35, 39].

In Fig. 2b, the system shows lasing at the band edge similarly as reported before in almost identical systems but pumping over the whole array, see [35, 36, 37] and references therein. Here the band edge is lower in energy than $\xi_{\mathrm{abs}}$, and thus absorption-emission cycles causing the red-shift are suppressed near the band edge. The linewidth of the emission peak is around 1.5 nm, which is narrower than the natural linewidth of the SLR modes (3.8 nm). This, together with a non-linear pump dependency with a clear threshold, is characteristic for lasing (see Extended Data Fig. 2).

The situation changes drastically when the band edge is matched with $\xi_{\mathrm{abs}}$, Fig. 2c. The emission at the band edge energy is negligible at the beginning of the array. Now the array periodicity is chosen so that absorption-emission cycles are possible until the band edge energy and thermalization may take place. A macroscopic occupation emerges when the thermalizating population reaches the band edge. In the following, we show that this phenomenon fulfills the main characteristics of a BEC: population following Bose–Einstein distribution with macroscopic occupation of the ground state and a thermalized tail, sharp increase in temporal coherence as well as build up of spatial coherence.

The difference between lasing and BEC is clearly shown in the evolution of the relative intensity at the band edge, and that of the rest of the population at higher energies, see Fig. 2f-g. In the lasing case, a lasing peak is generated at the pump spot (for spectra collected from the part of the array that is under the pump spot, see Extended Data Fig. 3) and then spreads through the array. In contrast, in the BEC case the population at the band edge is initially vanishing, and the thermalizing population is crucial for the emergence of the condensate peak. Thus we have a direct observation that the population accumulated to the ground state has emerged via thermalization. At around 150 $\mathrm{\mu m}$ the band-edge population becomes dominant, and after 200 $\mathrm{\mu m}$ the relative intensities of the band edge and the rest of the population stabilize. This is the area where we observe a room-temperature Bose–Einstein distribution. With the thermalization mechanism we use [22], quasi-equilibrium can be reached if multiple emission-absorption cycles take place [40, 41]. In our lasing case, absorption is suppressed, and therefore the thermal tail is vanishing (Fig. 2b,f and Extended Data Fig. 3). In the BEC case, the SLR excitations that form the condensate around 150 $\mu$m and thereafter are precisely those that have undergone several absorption-emission cycles (Fig. 2c,d,g,h). Therefore a Bose-Einstein (BE) distribution in our system has its origin in thermalization and condensation phenomena. This is important in order to distinguish from mere coexistence of a sharp peak and and a Boltzmann tail as observed for photonic lasing in semiconductor microcavities [42].

To further confirm that the observed luminescence originates from the SLR modes, we performed the measurements also in the momentum ($k$) space. As shown in Fig. 3a-b, the population indeed follows the SLR dispersion down in energy until the band edge. It follows one dispersion arm in Fig. 3a because, in the beginning of the array, only one SLR propagation direction is possible.

Transition as function of pump fluence Transition to BEC should occur in a non-linear manner at a critical density or temperature. Fig. 3c shows the measured population distributions for the BEC case at different pump fluencies. Using the relation between the critical density and the thermal de Broglie wavelength $\lambda_{T}=\sqrt{2\pi\hbar^{2}/(m_{eff}k_{B}T)}\simeq 10$ $\mathrm{\mu m}$ for a finite two-dimensional (2D) system, namely $n_{c}=2/\lambda_{T}^{2}\ln(L/\lambda_{T})$ [43, 9], we obtain critical density of the order $10^{10}/$m2 at $T=300$ K ($\hbar$ is the Planck’s constant and $k_{B}$ the Boltzmann constant). Finite systems with 2D linear or 1D parabolic dispersions give the same order of magnitude, and 1D linear an order of magnitude smaller. For the system size $L$, we used the measured coherence length of the SLR mode in the absence of molecules, $L=2\pi/\Delta k=45$ $\mathrm{\mu m}$. It was obtained from $\Delta k$ of the dispersion and confirmed by a double slit experiment, see Extended Data Fig. 4. Accurate determination of the density of photons at our sample is difficult, but the intensity of emitted photons gives it a lower bound. The lower bound is found to be higher, $10^{11}/$m2, than the critical density estimate (see Supplementary Information).

At small pump fluencies, the measured spectrum in Fig. 3c follows Maxwell–Boltzmann distribution (see Methods), but a peak emerges near the band edge when pump fluence is increased. After exceeding the threshold pump fluence, a spectrally sharp peak appears in the band edge with a long thermalized tail. The observed linewidth (1.5 nm) is narrower than the natural linewidth of the SLR modes, see also Extended Data Fig. 5. A decrease of the linewidth by around a factor of two is a direct indication of increase of the temporal coherence. Remarkably, the BE distribution fits excellently to the measured spectra in Fig. 3c, with $T=269\pm 67$ K, close to room temperature. Fits to the BE distribution for lower pump fluencies are shown in Extended Data Figure 6a. The distribution is an average over ten measurements, for details see Methods. In contrast, for the sample exhibiting lasing, the thermalized tail is vanishing and the fit to the BE distribution gives a temperature of only 25 K, see Extended Data Fig. 2.

Comparison to theory The interpretation of our experimental observations is supported by a rate-equation approach [40], see Methods. The computational spectral evolution in time-domain maps into the spatially measured spectra via group velocity of the SLR excitations. The experimentally observed evolution of the population, Fig. 2c and Fig. 2g, is in a good agreement with the computational result shown in Fig. 2d and Fig. 2h. For comparison between experiments in Fig. 2a-b and the theory, see Extended Data Fig. 7. The observed thermalization times are typically 0.5–1 ps. A pump-probe experiment further confirms the picosecond time-scale of the dynamics (Extended Data Fig. 8). The pump dependence of the population distribution (Fig. 3c-d) and the thermalization time dependence on concentration (Fig. 2e) show good agreement between the measurements and the model as well.

Off-diagonal long-range order

BEC is associated with spontaneous breaking of the U(1) symmetry via emergence of a well-defined global phase. In our case, this would correspond to a buildup of spatial coherence. The size, dispersion and dimension of the system affect BEC in profound ways. First, the dispersion and dimension determine the occupation probabilities of excited states and thereby the possibility of BEC [1, 43]. Second, fluctuations allow only quasi-long-range order in dimensions smaller than three [44]. It is thus essential to inspect the dimensionality, size, and dispersion of our system. SLR modes are known to be tightly confined in the third dimension, but whether our condensate is 1D or 2D is not determined a priori. By direct 2D $k$-space imaging we have confirmed that the condensate is confined in two dimensions, see Fig. 3e. Above we have shown in detail how the system thermalizes when the population is propagating along the $y$-axis (polarization along $x$), following the SLR dispersion in $k_{y}$. The dispersion in the $x$-direction is parabolic, and the dispersion has a 2D minimum at ${\bf k}=0$ ($\Gamma$-point) (Extended Data Fig. 9). The dispersions (an example shown in Fig. 1c) at high energies are linear, but around the band bottom, where the condensation occurs, they are flat and can be approximated for instance by a parabola. For 1D and 2D parabolic and 1D linear dispersions, a BEC is possible only in a finite system. Our system is finite: the effective size is given by the SLR mode coherence length (45 $\mathrm{\mu m}$). Although true long-range order is a subtle issue in low dimensions [44], especially in non-equilibrium systems [45], off-diagonal long-range order over a finite area may emerge at the critical density.

We measured the spatial coherence by a Michelson interferometer setup with one of the mirrors replaced by a retro-reflector that inverts the image along the y-axis, see Methods. Combining the real-space image of the sample with its mirror image reveals the phase correlations between different parts of the sample. From the measured pattern shown in Fig. 4a, we observe that the condensate is spatially coherent over distances of at least 90 $\mathrm{\mu}$m, confirmed by Fig. 4c. This is twice as long as the natural coherence length of the SLR modes (45 $\mathrm{\mu}$m), demonstrating that the degree of spatial coherence has increased. With increasing pump fluence above the threshold value, we see a gradual build-up of the degree of spatial coherence as shown in Fig. 4b.

Conclusions

We studied nanoparticle arrays with dye molecules acting as a thermal bath, and monitored propagation of the surface lattice resonance (SLR) excitations after pumping. Fit to the BE distribution at room temperature, spectral narrowing, non-linear behaviour at critical density, as well as the onset and growth of spatial coherence were observed in our finite-size, open-dissipative system. Difference to lasing was demonstrated by suppressing the possibility of the population reaching the band edge via thermalization. The evolution of the band edge peak and thermalizing population revealed that lasing occurs at the pump spot while the BE distribution emerges from thermalizing population. These observations constitute evidence for a BEC of surface plasmon polariton lattice modes.

The sub-picosecond dynamics in our system is faster than the typical condensation dynamics in magnon condensates [20] (100 ns), photon condensates[22, 29, 23] (10 ps – 1 ns) and most inorganic polariton condensates [9, 10, 12] (1–100 ps). It is comparable to those in the fastest inorganic and organic polariton systems [16, 18, 19] where sub-picosecond decay and coherence times have been observed. Our experiment allows monitoring the actual evolution of the thermalization and condensate formation in detail, as exemplified by Fig. 2, something that has earlier only been possible in longer time-scales.

Strong coupling has been achieved in plasmonic systems at room temperature [38, 46], and the SLR modes can hybridize with the molecules [47]. This may provide effective interactions and thus a new arena for photonic quantum fluids. The lattice geometry can be changed into more complex ones, such as hexagonal [48], in search for degeneracy (Dirac) points. Customizable lattice geometries are exciting, especially because recent discoveries of topological phenomena [49] have accentuated the significance of the band structure. Symmetries can be broken, for instance, by nanoparticle shape or use of magnetic material [50] to explore how condensation may be transformed in a topologically non-trivial system. Finally the new type of BEC we have created is highly accessible technologically: it is based on on-chip systems that can be scaled up by using nanolithography and microfluidistics, and operate at room temperature.

Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper.

Methods

Sample fabrication

Square arrays of rod shaped gold nanoparticles were fabricated on a borosilicate substrate by electron-beam lithography and subsequent electron-beam evaporation of 2 nm of titanium and 50 nm of gold. The array periodicities ranged from $p_{x}=p_{y}=$ 580 nm to 610 nm. The width of the particles was fixed to 100 nm and the length was chosen to be 65 $\%$ of the period $p_{y}$ to provide a sufficient band gap in the dispersion. The IR-792 molecules were dissolved into a solution of 1:2 (dimethyl sulfoxide):(benzyl alcohol). The time-resolved photoluminescence (TRPL) of the molecule does not show any wavelength dependence (Extended Data Fig. 11) that would complicate the analysis of the results.

Measurement setup

The $x$-axis of the 2D charge coupled device (CCD) camera of the spectrometer was used to resolve the wavelength. Depending on the measurement scheme, the $y$-axis was used to resolve either 1) the angle of the transmitted/emitted light or 2) the $y$-position along the array. The spectra were collected over 300 ms integration time, that is, averaged over 300 laser pulses.

Angle-resolved measurements

Angle-resolved transmission measurements (as in Fig. 1c) were obtained by focusing light from a halogen lamp onto the sample and collecting the transmitted light with an objective (10X, 0.3 NA). The back focal plane of the objective was focused to the entrance slit of the spectrometer, and thus each position along the $y$-axis of the CCD corresponded to a particular angle of transmitted light. The $E(k)$ dispersions were subsequently calculated from the angle and wavelength-resolved spectra as $E=hc/\lambda_{0}$ and $k=k_{0}sin(\theta)$. Here $k_{0}=2\pi/\lambda_{0}$ and $\theta$ is the angle with respect to sample normal. The dispersions show a band edge, and we have confirmed that the band-edge energy remains the same (with accuracy of four digits in eV) throughout the array (Extended Data Fig. 12).

Spatially specific k-space images of the sample emission (such as in Fig. 3), were measured in a similar fashion, but the sample was pumped with a pulsed femtosecond laser (100 fs pulse duration, 750 nm center wavelength, 1 kHz repetition rate, $\mathrm{\mu}$J pulse energy), with a flat pulse profile. Pulsed excitation allows for high photon densities necessary for quantum degeneracy, while simultaneously avoiding bleaching of the organic molecules. Further, an iris was placed on the image plane of the sample to spatially filter only the light emerging from a specific region of the sample. For the spatial coherence measurements, the iris was replaced with a double slit having a 90 $\mathrm{\mu}$m inter-slit distance and 30 $\mathrm{\mu}$m slit width. For a schematic, see Extended Data Fig. 13a.

Spatially resolved measurements

For spatially resolved spectra such as in Fig. 2, the real space image of the sample was focused to the spectrometer entrance slit. Thus the $y$-axis of CCD corresponded to a particular position along the $y$-axis of the array while the $x$-axis corresponded to the wavelength. The magnification of the real space image was 4x at the entrance slit and the entrance slit width was 300 $\mathrm{\mu m}$ so that we collected the luminescence from 3/4 of the 100 $\mathrm{\mu m}$ array (in $x$-direction). For a schematic, see Extended Data Fig. 13b.

Michelson interferometer

Spatial coherence properties were studied with a modified Michelson interferometer configuration, for a schematic see Extended Data Fig. 13c. The luminescence of the sample was first passed through a 900 nm (1.378 eV) long pass filter and then was directed to a non-polarizing beam splitter. One output arm of the beam splitter consisted of a hollow roof L-shaped mirror used to invert y to –y upon reflection. The other arm consisted of a regular mirror. The two reflected beams were directed through the beam splitter onto a lens and then onto a CCD. The first-order correlation function describing the degree of spatial coherence reads $g^{(1)}(\boldsymbol{y},\boldsymbol{y}^{\prime})=\frac{\langle E^{*}(\boldsymbol{y})E(\boldsymbol{y}^{\prime})\rangle}{\langle E^{*}(\boldsymbol{y})\rangle\langle\langle E(\boldsymbol{y}^{\prime})\rangle}$ (here $E^{*}(\boldsymbol{y})$ is the electric field at point $\boldsymbol{y}$). The contrast is obtained directly from the interference pattern as

[TABLE]

Here, $I(\boldsymbol{y})$ is the light intensity and $I_{\mathrm{max}}$ and $I_{\mathrm{min}}$ are the maximum and minimum intensities of the fringes, respectively. Note that in our case the contrast is expected to decrease for increasing $y$ due to the decay ($I(\boldsymbol{y})<I(-\boldsymbol{y})$) of the SLR population, therefore the contrast gives only a lower limit to the value of $g^{(1)}(\boldsymbol{y},-\boldsymbol{y})$.

Thermalization rate measurement

The dependence of thermalization rate on concentration (Fig. 2e) was determined by repeating the thermalization experiment of Fig. 2a with different concentrations (25, 50, 75, 100, and 150 mMol). From the intensity spectra (as Fig. 2a-c), we then obtained the slope of the red shift by fitting a line to the maxima of the intensity at each $y$ position. The slope of the red shift interprets as the rate of the thermalization process by scaling the distance with the SLR-mode group velocity. The slopes were measured from the energy regime that is far from the band edge, i.e. from the linear part of dispersion. Various periodicities were used in the experiments. The data in Fig. 2e shows the mean value of several (3–7) measurements for each concentration, with error bars showing the standard deviation within each data set. Note that the linear fits for thermalization rate are valid in a limited energy range. This is evident from Extended Data Fig. 14 where the thermalization and emergence of BEC is shown for a period $p=600$ nm sample with a broader energy range (at the expense of spectral resolution).

Theoretical model

We used a rate-equation model [40] that describes population dynamics but neglects correlations relevant for strong coupling. In the model, an ensemble of spectrally broadened two-level dye molecules is pumped by a laser pulse. The dye molecules interact with discretized, lossy SLR modes by emission and re-absorption. Measured emission and absorption profile taken from literature, as well as realistic parameters for the loss rates and lifetimes were used in the model. Parameters in Eq. 2 of Supplemental Information, used in the simulations are the following: SLR decay rate $\kappa=1.08\times 10^{13}$ 1/s, number of molecules is varied $N=1.18...7.08\times 10^{7}$ ($N=2.34\times 10^{7}$ corresponds to experiments with concentration of 50 mMol), coupling constant $g=10^{4.5}$, pump amplitude is varied $P=5.01\times 10^{10}...13.80\times 10^{10}$, and spontaneous decay rate for the molecules $\gamma_{m}=10^{9}$ 1/s. With these parameters, the fraction of molecules that is excited by the pump is small, therefore there is a large number of molecules in thermal equilibrium available for re-absorption processes. The main approximations compared to a full microscopical description are that a phenomenological parameter (the coupling $g$) is used for the overall strength of absorption and emission, and the field profiles are assumed uniform. There are no computational states in the region of the energy band gap of the dispersion, which produces an empty area in Fig. 2d and Fig. 3d. The calculation of the thermalization time as a function of molecular concentration, Fig. 1e, constitutes a self-consistency check for the correspondence between the experiment and the model. Namely, we varied the concentration while keeping all other parameters fixed, and obtained the good agreement between the model and the measurements shown in Fig. 2e. For detailed formulation of the model, see Supplementary Information.

Fits to Bose–Einstein and Maxwell–Boltzmann distributions

For the measured population distributions, examples shown in Fig. 3 and Extended Data Fig. 6, fits to Bose–Einstein (BE) and to Maxwell–Boltzmann (MB) distributions were made using non-linear least squares fitting. In a few meV vicinity of the $\Gamma$-point the dispersion is quite flat both in $x$- and $y$-directions, and for higher energies it is essentially linear but anisotropic in the two dimensions as shown in Extended Data Fig. 9. The density of states used in fitting was obtained numerically from such dispersion. Above threshold, excellent fits to the BE distribution around room temperature were obtained while the MB distribution did not fit well. Below threshold, the MB distribution gives good results at $T=300$ K. To obtain statistics for the temperature of the condensate, we made the BE distribution fit to data obtained by averaging over ten separate measurements. All the measured population distributions were taken using the corresponding pump fluencies and from the same region of the array (the red circle in Fig. 3). The resulting fit to the BE distribution gives temperature of $T=269\pm 67$ K, error limits showing 95% confidence bounds for the fit. With 99% confidence bounds the result is $T=269\pm 89$ K. In contrast, for the sample exhibiting lasing the fit to the Bose–Einstein distribution gives a temperature of only 25 K, see Extended Data Fig. 2, reflecting the non-equilibrium character of lasing and lack of thermalization to the environment.

Acknowledgements We thank Miikka Heikkinen, Dong-Hee Kim, Robert Moerland and Marek Nečada for useful discussions. This work is dedicated in memory of D. Jin and her inspiring example. This work was supported by the Academy of Finland through its Centres of Excellence Programme (2012–-2017) and under project NOs. 284621, 303351 and 307419, and by the European Research Council (ERC-2013-AdG-340748-CODE). This article is based upon work from COST Action MP1403 Nanoscale Quantum Optics, supported by COST (European Cooperation in Science and Technology). K.S.D. acknowledges financial support by a Marie Skłodowska-Curie Action (H2020-MSCA-IF-2016, project id 745115). Part of the research was performed at the Micronova Nanofabrication Centre, supported by Aalto University. Triton cluster at Aalto University was used for the computations.

Author contributions P.T. initiated and supervised the project. T.K.H., A.J.M., R.G., and A.I.V. performed the experiments. A.J.M., T.K.H., and A.I.V. analysed the data. T.K.H., K.S.D., and H.T.R. built the experimental setup. A.J.M., J.-P.M., and A.J. performed the theoretical modeling. R.G. fabricated the samples. All authors discussed the results. P.T., A.J.M., and T.K.H. wrote the manuscript together with all authors.

Additional information Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints. Correspondence and requests for materials should be addressed to P.T.

Competing financial interests The authors declare no competing financial interests.

Supplementary Information

Rate equations

Our theoretical rate equation model is based on the approach introduced in [1, 2] where a microscopic quantum model was developed to describe the dynamics of photon condensates in cavities. Our system consists of discretized SLR modes which are coupled with $N$ dye molecules. The model assumes that each molecule is a two-level system coupled to a harmonic oscillator that describes the rovibrational degrees of freedom within a molecule. For such a system the microscopic Hamiltonian is of the form

[TABLE]

Here $a_{\textbf{k}}^{\dagger}$ is the bosonic creation operator of the SLR mode of momentum k, $\sigma_{i}^{z}$ and $\sigma_{i}^{\pm}$ are the Pauli operators describing the two-level structure of the $i$th molecule and $b_{i}^{\dagger}$ is the bosonic creation operator corresponding to the vibrational ladder of the $i$th molecule. Furthermore, $\omega_{\textbf{k}}$ is the energy of the SLR mode of momentum k, $\Omega$ is the vibrational energy of the molecules, $g_{i}$ describes the coupling strength between the $i$th molecule and the SLR modes, $\omega_{m}$ is the energy of the two-level systems and $X$ describes the spectral broadening of the molecules due to the vibrational degrees of freedom. In Eq. (2) we have assumed that electric field profiles of SLR modes are uniform which is an approximation. In plasmonic lattices the electric fields can have rather complicated profiles and quantizing such fields is extremely challenging. For details of deriving the above Hamiltonian, see for instance [3, 4, 5, 6, 7].

Our system is a highly open quantum system as the SLR modes are extremely lossy, the molecules are externally pumped by a laser beam, and the rovibrational degrees of freedom have small relaxation times. We treat our system by writing a master equation for the time derivative of the reduced density matrix for the SLR modes, $\frac{d\rho}{dt}$, and deploying the standard Born-Markovian approximation for the pump, SLR losses and rovibrational states [1, 2, 3]. With $\frac{d\rho}{dt}$ we can solve the rate of change of the SLR populations $n_{\textbf{k}}=\langle a_{\textbf{k}}^{\dagger}a_{\textbf{k}}\rangle$ by writing $\frac{d}{dt}\langle a_{\textbf{k}}^{\dagger}a_{\textbf{k}}\rangle=\textrm{Tr}[a_{\textbf{k}}^{\dagger}a_{\textbf{k}}\frac{d\rho}{dt}]$ and the fraction of excited molecules $p_{e}=\langle\sigma_{i}^{+}\sigma_{i}^{-}\rangle$ by having $\frac{d}{dt}\langle\sigma_{i}^{+}\sigma_{i}^{-}\rangle=\textrm{Tr}[\sigma_{i}^{+}\sigma_{i}^{-}\frac{d\rho}{dt}]$. The rate equations used in this work become:

[TABLE]

Here $\kappa(\textbf{k})$ is the decay rate of the SLR mode of momentum k, $\Gamma^{\pm}(\textbf{k})$ are the molecular emission ($+$) and the absorption ($-$) coefficients for a given SLR mode with momentum k, $P(t)$ is time-dependent pump power and $\gamma_{m}$ describes the spontaneous non-radiative decay of the molecular two-level systems from the upper energy level to the ground state. One can derive microscopic expressions for $\Gamma^{\pm}(\textbf{k})$ but in this work we used real emission and absorption profiles of the molecule.

The derivation of Eqs. (Rate equations) requires multiple approximations. First, we assumed that the correlator term $\langle a_{\textbf{k}}^{\dagger}a_{\textbf{k}}\sigma_{i}^{+}\sigma_{i}^{-}\rangle$ can be approximated as $\langle a_{\textbf{k}}^{\dagger}a_{\textbf{k}}\sigma_{i}^{+}\sigma_{i}^{-}\rangle\approx\langle a_{\textbf{k}}^{\dagger}a_{\textbf{k}}\rangle\langle\sigma_{i}^{+}\sigma_{i}^{-}\rangle$. This is likely a valid approximation as we are considering the weak-coupling regime. The approximation makes it possible to describe the system completely with Eqs. (Rate equations) without a need to write extra rate equations for $\langle a_{\textbf{k}}^{\dagger}a_{\textbf{k}}\sigma_{i}^{+}\sigma_{i}^{-}\rangle$ and other higher-order correlators.

Furthermore, we assumed that all the molecules behave similarly. This can be justified by recalling that all the molecules interacting with plasmonic modes are in the proximity of the nanoparticles and thus the volume enclosing the active molecules is small. Therefore we assume that the population inversion $p_{e}$ and the coupling coefficient $g_{i}=g$ are the same for all molecules which means we need only one rate equation to describe the dynamics of all of them.

Spatial coherence measurement

To study the emergence of spatial coherence, we imaged the condensate by a Michelson interferometer with one of its arms replaced by a hollow roof L-shaped mirror which inverts the image along the $y$-axis. In this configuration we measure the phase coherence between points at $y$ and $-y$. The intensity of the interference pattern between two plane waves is given by [8]

[TABLE]

where $A$ is a constant related to the speed of light and $\boldsymbol{E}_{1}=E_{1}\textrm{exp}(i[\boldsymbol{k}_{1}\boldsymbol{r}-\omega t+\phi_{1}])$ and $\boldsymbol{E}_{2}=E_{2}\textrm{exp}(i[\boldsymbol{k}_{2}\boldsymbol{r}-\omega t+\phi_{2}])$ are the electric fields from the arms of the Michelson interferometer and $\boldsymbol{E}_{1}\boldsymbol{E}_{2}\cos{\theta}$ is the interference term responsible for the fringe spacing, orientation and initial phase. The first-order correlation function, describing the degree of spatial coherence, between the points $\boldsymbol{r}$ and $\boldsymbol{r}^{\prime}$ is

[TABLE]

When the mirror image is taken centrosymmetrically (in our case $\boldsymbol{y}^{\prime}=-\boldsymbol{y}$), the correlation function is obtained directly from the visibility of the interference fringes, as

[TABLE]

Here, $I(\boldsymbol{r})$ is the light intensity and $I_{\mathrm{max}}$ and $I_{\mathrm{min}}$ are the maximum and minimum intensities of the fringes, respectively. Note that in our case the contrast is expected to decrease for increasing $y$ due to the decay ($I(\boldsymbol{y})<I(-\boldsymbol{y})$) of the SLR population.

Estimating photon densities

Critical photon density for BEC

We calculate the critical photon density for BEC in a finite 2D system [9, 10, 11, 12]. Assuming the SLR mode dispersion to be parabolic around the band bottom,

[TABLE]

we obtain an estimate for the effective mass $m_{eff}$. Here, $\hbar$ is the reduced Planck’s constant $h/2\pi$ and $k$ the momentum. Fitting Eq. (7) to the band bottom gives $m_{eff}=1.7...$ $2.6\times 10^{-37}$ kg. The values $1.7\times 10^{-37}$ kg and $2.6\times 10^{-37}$ kg are obtained by two equally reasonable fits where different range of wavevectors around k = 0 are selected as the band bottom. With the standard 95 % confidence level, the lower and upper error limits (standard 95 % confidence level) for these values are $1.3\times 10^{-37}$ and $4.3\times 10^{-37}$, respectively. The respective 99 % confidence bounds for the values are $1.2\times 10^{-37}$ and $5.8\times 10^{-37}$. Critical density ($n_{c}=N/V$) is given by

[TABLE]

where $L$ is the size of the 2D box. We take the measured coherence length of the SLR modes in absence of molecules ($45\mathrm{\mu m}$) as the box width. Thermal de Broglie wavelength is defined as

[TABLE]

where $k_{B}$ is the Boltzmann constant. Now, using Eq. (8) and the effective mass, the critical photon density is between $n_{c}=3.0...5.3\times 10^{10}$ 1/m2. Within 95 % confidence bounds for the effective mass, the critical photon density is between $n_{c}=2.1...10.0\times 10^{10}$ 1/m2. With the 99 % confidence bounds, the critical photon density is between $n_{c}=1.8...14.6\times 10^{10}$ 1/m2. Repeating such calculation for 2D linear and 1D parabolic dispersions gives the same order of magnitude, $n_{c}\simeq 10^{10}$ 1/m2, while 1D linear gives $n_{c}\simeq 10^{9}$ 1/m2 (finite size is required in the 1D linear/parabolic cases but not in 2D linear one). Therefore the estimate $n_{c}\simeq 10^{10}$ 1/m2 is reasonable for our real dispersion which is a combination of linear and flatter regions, and is anisotropic. One has to also bear in mind that these are equilibrium values while our system is open-dissipative and can at best reach a quasi-equilibrium.

Lower bound for the photon density in the experiment

As described in Methods, the sample luminescence was measured with the 2D charge coupled device (CCD) camera. The CCD intensity reading was calibrated against a 633 nm HeNe laser, taking into account the wavelength-dependent quantum yield of the CCD. The emission was collected through the spectrometer entrance slit whose width was 3/4 of the array width on the sample.

By summing the collected intensity over the part of the sample where the condensate is observed, we get the total number of counts in the CCD to be around $4.98\times 10^{4}$. Taking into account the overall efficiency of the CCD (limited by instrument response, quantum yield, etc.), we measured that each CCD count requires 20 photons. The luminescence is measured only from one side of the sample, and we assume the sample to radiate equally to both sides, so the actual photon number in the system is twice the observed. Therefore, the lower limit for the total number of photons in the BEC is around $1.99\times 10^{6}$.

As the integration time in the luminescence measurements was 300 ms, and the laser pulse frequency is 1 kHz, we can deduce the number of photons in the BEC per pulse to be $6.63\times 10^{3}$ . Dividing this by the area from where the emission was collected, 75 $\times$ 100 $(\mathrm{\mu m})^{2}$, the lower bound for the number density of photons in the BEC is around $8.8\times 10^{11}$ 1/m2.