Supplemental Material for: Quantum fluids of light in all-optical scatterer lattices

S. Alyatkin,1, 2 H. Sigurdsson,1, 2, 3 A. Askitopoulos,1, 2 J. D. Töpfer,1, 2 and P. G. Lagoudakis1, 2, 3, ∗ 1Center for Photonics and Quantum Materials, Skolkovo Institute of Science and Technology, Moscow, Russia 2Laboratories for Hybrid Photonics, Skolkovo Institute of Science and Technology, Moscow, Russia 3School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, United Kingdom (Dated: August 24, 2021)

In this section we provide additional data on the measured polariton photoluminescence (PL) in the diffractive scatterer regime using pump spots arranged in a square lattice with lattice constant D = 22.5 µm (see Fig. S1) and Lieb lattice D = 22.3 µm (see Fig. S2). We also present complementary data on the inverse Lieb lattice in Fig. 4 in the main manuscript (see Fig. S3). The full-width-half-maximum (FWHM) of the Gaussian laser spots exciting the condensates is ≈ 2 µm. Just like in the main text, we define the real space horizontal and vertical coordinates as r = (x, y), and in momentum space k = (k x , k y ).
Figures S1a and S1b show the PL in real space and energy-momentum space (along k x = 0) approximately at threshold (P ≈ P thr ) for a square arranged lattice of widely spaced Gaussian pump spots. At this power no significant phase locking occurs between the condensates as can be evidenced from the absence of interference fringes in the spatial PL and the broad energy distribution of polaritons in the lattice dispersion. In Figs. S1c and S1d we show the same measurements above condensation threshold (P ≈ 1.3P thr ) which result now in clear appearance of interference fringes between the radiating condensates and a greatly narrowed energy distribution peaking around 1.4478 eV. We note that Figs. S1b and S1d are plotted on a logarithmic scale. We observe multiple illuminated high energy bands in the polariton dispersion verifying that polaritons in the system are experiencing crystal scattering in the optical lattice.
In Fig. S2 we provide similar measurements to those in Fig. 2a-e in the main manuscript but using a larger lattice Figure S1. Optical lattice of 25 polariton condensates arranged in a square geometry with a lattice constant set to D = 22.5 µm. (a,c) Real space and (b,d) energy-resolved momentum space polariton PL at (P ≈ P thr ) and above (P ≈ 1.3P thr ) threshold respectively. Figure S2. Optical lattice of 96 polariton condensates arranged in a Lieb geometry (a) at threshold with a lattice constant set to D = 22.3 µm. The inset in (a) schematically shows sublattices denoted with letters A,B,C forming the Lieb lattice. (b) Real space polariton PL for the lattice pumped at P ≈ 1.1P thr . Driving the system above threshold leads to the condensation on A and C sublattices, see inset shown in linear colorscale. (c) Momentum space polariton PL for P ≈ 1.1P thr . Semi-transparent thin lines denoted with yellow letters "d" and "e" correspond to kx = 0 and kx = 2π/D, for which the energy-resolved momentum space PL shown in (d),(e) respectively.
constant of D = 22.3 µm. The results presented in Fig. S2 show that by changing the lattice constant in the refractive (scatterer) regime we can change the qualitative features of the polariton PL above threshold. Here we observe destructive interference inside the cells of the lattice (i.e., where pump spots are absent) whereas in Fig. 2 we saw constructive interference. This result can be intuitively understood from the fact that polaritons are condensing into a favorable Bloch mode of the optical lattice by maximising their mutual interference (overlap) over the pump gain region. Changing the lattice constant D naturally affect the interference (i.e., more wavelengths fit between pump spots) and one observes distance-periodic changes in the polariton PL pattern as different Bloch modes move in and out of the lattice gain bandwidth.
In Fig. S3 we provide supplemental experimental results for the inverse (conventional tight binding) Lieb lattice hosting P -flatbands corresponding to Figs. 4e-g in the main text. Figure S3a shows the energy-resolved reciprocal space PL taken at k x = 0 where yellow arrows indicate the location of P -and S-bands. Figure S3b is measured at k x = 2π/D where the two horizontally dashed lines, denoted with yellow letters "(c)" and "(d)" correspond to the momentum space isoenergy planes PL shown in Figs. S3c and S3d. We note that in each panel the PL is independently normalised. The data shown in Fig. S3c is the same as shown in Fig. 4l in the main manuscript and corresponds to polariton condensation into the M symmetry point of the P -flatband. Unfortunately, the limited energy resolution of our experimental setup and the small non-stimulated population of polaritons in surrounding band structures prevents us from accurately resolving the P -flatbands in energy-momentum space. We additionally present the momentum space PL corresponding to the lowest P -band, forming an peculiar square shape pattern along the Γ ↔ X edges in the outer Brillouin zones. The lowest P -band is completely flat along the Γ ↔ X edges but possesses strong dispersion along the X ↔ M ↔ Γ edges as shown in Fig. 4i in the main manuscript. The same type of PL is also observed in our calculations (see Fig. S7e).

S2. THEORY OF POLARITON DYNAMICS
We consider the two-dimensional (2D) Schrödinger equation which describes a polariton in the static potential landscape designed by our optical non-resonant lasers.
Here,Ê is the dispersion of the lower polariton branch. The upper polariton branch is neglected because it is both far away in energy and weakly populated by our experimental excitation scheme. Given the much larger effective mass of the excitons compared to the cavity photons we can set their dispersion to zero and the lower polariton dispersion in momentum space becomes 1 ,Ê whereÊ φ (k) is the cavity photon dispersion and Ω = 12 meV is the Rabi energy of our sample 2 . The longitudinally confined photons obey the dispersion relation, where λ = 860 nm is the cavity resonance, n = 3.556 is the cavity refractive index, ∆ = −5 meV is the cavity photon detuning from the exciton resonance, and γ = 1/5.5 ps −1 is the cavity photon lifetime. The real part of the calculated dispersion is shown with a white dashed line plotted against the experimentally retrieved dispersion in Fig. S4a. In Fig. S4b we show a schematic of the Lieb lattice and the three square sublattices denoted by their label A, B, C.
The potential V (r) is directly proportional to the incident non-resonant laser intensity and can be written as a superposition of spatially separated Gaussians, Here, N is the number of pumps in question and x n , y n are their coordinates, V 0 denotes the complex potential strength. Re(V 0 ) > 0 and Im(V 0 ) > 0 denote the blueshift and particle gain experienced by polaritons respectively due to their interactions with the laser generated background of non-condensed particles, and w = 1.274 µm is the RMS width of the pump spots corresponding to 3 µm FWHM. The width is slightly larger than that of the lasers (2 µm FWHM) to account for the finite diffusion of excitons in the reservoir from the laser spot. The dispersion of Eq. (S1) can be numerically obtained through Monte-Carlo methods where an average is performed over the time evolution of many random initial realizations of the particle wavefunction Ψ(r, t = 0). This method is particularly useful for potentials which lack the usual symmetries required by Bloch's theorem to work. It permits calculation over which frequency and momentum componentsψ(ω, k) survive the longest in the system (i.e., decay the slowest) in presence of optical gain at the pump spots [i.e., Im(V 0 ) > 0]. Results of this method are shown in Figs. S5-S7 where we extract numerically the dispersion from the dynamics of many random initial realizations ofψ(ω, k) in Eq. (S1) in several lattice geometries relevant to the study. Overall we find good agreement between experimental observations and our calculations across all lattice geometries.

S3. BLOCH THEOREM
In this section we investigate the dispersion belonging to Eq. (S1) by applying Bloch's theorem. The following equations were used to calculate the bands plotted in Figs. 2i and 4i in the main text. For simplicity, we apply a parabolic approximation by writing,Ê where the effective mass of the lower polaritons is m = 0.32 meV ps 2 µm −2 found by fitting the dispersion in Fig. S4a. The potential is assumed to be infinite and periodic such that V (r) = V (r + D) where D = n 1 a 1 + n 2 a 2 is the translational symmetry vector defined in the bases of primitive lattice vectors a 1,2 for some integers n 1,2 . We can then apply Bloch's theorem where we write the polariton single particle wavefunction in the factorised form of crystal momentum q = (q x , q y ) and Bloch states in the nth band u n,q (r), Here, u n,q (r) = u n,q (r + D). Substituting Eq. (S6) into Eq. (S1) the Bloch eigenvalue problem then reads,    We can diagonalise Eq. (S7) by expressing the Gaussian potential as a Fourier series in the basis of the reciprocal lattice vectors defined through G · D = 2πn for n ∈ Z, (S8) Figure S8a shows a portion of the calculated band structure for a Lieb arrangement of the Gaussian spots (see Fig. S8f) with a lattice constant of D = 22.3 µm with some features of interest marked with magenta colored boxes. Figure S8b shows the real space intensity of the Bloch state |u nc,q=0 | 2 , marked by red circle corresponding to the experimentally observed pattern shown in Fig. S2. A portion of the Gaussian pump locations are marked by red squares in both Figs. S8b and S8f for clarity. In Fig. S8a we observe regions in the band structure with almost linear-and flat-like bands (magenta squares zoomed in Figs. S8c-e). Although such exotic points were not experimentally investigated in this study, their theoretical prediction enhances the versatility and importance of the investigated scatterer Lieb lattice. Even exciting polaritons into such points via non-resonant methods is infeasible they can still be excited using resonant lasers. This can potentially open up investigation into nonlinear transport properties of nearly massless or slow-light polaritons. We also observe that the calculated Bloch state shown in Fig. S8b shows complete destructive interference at sublattice B whereas in the experimentally measured spatial PL in Fig. S2b this destructive interference is not complete. This could mean the the condensate is populating more than one Bloch state of the lattice since Fig. S8a clearly shows that multiple bands exist within a narrow energy range (i.e., many bands can be within the gain bandwidth of the lattice). The time-integrated measurements of our experiment would then instead show a weighted average over these multiple Bloch states which explains the discrepancy between linear Bloch's theory and our measurements. In order to get better agreement between measurements and theory, we introduce in the next section a generalised Gross-Pitaevskii model which describes the polarion dynamics above condensation threshold (i.e., in the nonlinear regime).

S4. GENERALISED GROSS-PITAEVSKII SIMULATIONS
We simulate the dynamics of the polariton condensate in the parabolic regime of the lower polariton branch using the two-dimensional driven-dissipative (i.e., generalised) Gross-Pitaevskii model 3 written for the condensate wavefunction, Ψ(r, t), and a rate equation for the exciton reservoir supplying particles into the condensate, n(r, t): Here, m is the polariton effective mass, α is the polariton-polariton interaction strength, R is the scattering rate of reservoir excitons into the condensate, γ is the rate of polariton losses through the cavity mirrors, Γ R is the decay rate of the reservoir excitons, g and G are the interaction strengths of polaritons with the exciton reservoir feeding the condensate (i.e., so-called bottleneck or active excitons), and high momentum photoexcited excitons P (r) (proportional to the laser profile and power) that do not scatter into the condensate, respectively. The final term R(Ψ) is an energy relaxation term 4,5 that assists ground state condensation and is taken proportional to the background density of bottleneck excitons, Here, λ denotes the energy relaxation efficiency, ξ is a fitting parameter, and µ(r, t) is the particle conserving local effective chemical potential of the condensate.
To take into account to finite diffusion of the reservoir excitons n(r, t) away from the laser pump spots we approximate the reservoir driving profile with a Lorentzian, (S12) Here, N is the number of pumps in question and x n , y n are their coordinates, 2L = 2 µm denotes the full width at half maximum of the reservoir driving field and the value P 0 denotes the laser power density. We point out that the band features previously analysed in Sec. S2 using lattices made up of arrranged Gaussian spots are not qualitatively changed here. While our model aims at being as simple as possible to produce the observed effects, more accurate simulations could possibly be achieved by introducing explicitly an exciton diffusion term into Eq. (S10) 6 and stochastic noise 7 . The lower polariton mass and lifetime are taken corresponding to the properties of our cavity: m = 0.32 meV ps 2 µm −2 , and γ = 1/5.5 ps −1 . The remaining parameters are taken similar to those used in previous works: α = 3.3 µeV µm 2 , g = α, G = 2g/Γ R , R = 0.8α, Γ R = 2γ, ξ = 1/4, and λ ξ = 0.004 µm 2 .
In each simulation we apply stochastic initial conditions to Ψ(r, t = 0) and run the simulation much longer than the characteristic timescales of the condensate dynamics. This typically corresponds to simulation times longer than 10 ns. The results of such simulations are presented in the insets of Figs. 3a-c and 3e in the main text where we have plotted the time-integrated (averaged) condensate density |Ψ(r, t)| 2 at the same lattice constants as in experiment for pump powers P 0 chosen to give the best match. The single exception is the data presented in the inset in Fig. 3d which was not obtained from Gross-Pitaevskii simulations. This was due to the lack of the simulated condensate stabilising in a state similar to the one observed in experiment for the chosen parameters. It is important to stress that this does not preclude such a condensed state existing in stable form in simulation for other parameters. However, we are able to find a good match using Bloch's theorem [see Eq. (S7)] where we plot in Fig. 3d the state |Ψ(r)| 2 = |u 3,q X (r)| 2 +|u 3,q X (r)| 2 with q X = (0, π/D) and q X = (π/D, 0) corresponding to the edges of the Brillouin zone.
In Fig. 3m in the main manuscript we qualitatively produce the experimental observation of gain guided condensates collapsing into trapped condensates at high power and small lattice constant through Gross-Pitaevskii simulations. Some discrepancy can be observed between simulation and experiment at low powers and small lattice constants where the simulation overestimates the gain guided polaritons. This discrepancy can appear from the lack of stochastic treatment (i.e., Wiener noise) in our equations which would smear out the simulated condensate PL at low powers close to threshold.

S5. FINITE P -FLATBAND CURVATURE
As shown in Fig. 4i and 4j in the main manuscript, there exist residual finite band curvature in the P -flatbands of the inverse optical Lieb lattice. As mentioned in the text of of main manuscript the dispersion (curvature) of the flatband appears around the M point and is ≈ 5 µeV different in energy from both the Γ or the X points (see  The source of this curvature is attributed to the finite potential depth and the Gaussian shape of the excitation beams (see Fig. S9a) which leads to deviation from the perfectly dispersionless bands predicted by the tight binding theorem. By increasing the depth of the potential minima (or conversely, the height of the potential maxima) in the lattice the flatband curvature is reduced as can be seen in Fig. S9c. The flatbands in Fig. S9c are also blueshifted (with respect to the other P bands) which comes from a slightly increased energy splitting between the P x and P y orbitals at sublattices A and C due to the Gaussian construction of our lattice.

S6. TIGHT BINDING TREATMENT
Here, we show qualitatively why condensation preferably occurs into a flatband P-orbital based by considering a tight-binding treatment of the P -orbital states. We consider only nearest neighbor coupling and negligible mixing with neighboring S-and D-orbitals. Being a scalar problem, we also do not consider the polarization degree of freedom of the polaritons. The bulk Hamiltonian in the basis of single particle, lattice site, creation and annihilation operators can be written,Ĥ = n,m n,m−1 )] + h.c.
Performing the standard Fourier transform of the single particle operators one can transform the problem from the (S14) Setting all diagonal terms to zero corresponds to all sites being in resonance and with equal on-site gain. In Fig. S10a we plot the dispersion from Eq. (S14) for J 1 = 0.4 + i0.04 and J 2 = J 1 /3. The blue-red colorscale indicates energies with low and high particle gain (imaginary part of the energy), respectively. In this case it shows that polaritons will survive the longest around the Γ point in the upper bands.
On the other hand, if the x and y P -orbitals become split in both real energy and gain at their respective sublattice sites then the dispersion dramatically changes. In Fig. S10b we show the calculated energies again for J 1 = 0.4 + i0.04 and J 2 = J 1 /3 but with x and y P -orbitals split at sublattices A and C. This is modeled by setting a,x = c,y = − 0 + iκ, a,y = c,x = 0 − iκ, and b,x = b,y = 0 where 0 , κ > 0. Such splitting arises from the slightly elliptically shaped A and C sublattice sites [see Fig. S7a]. Physically, κ captures the increased gain of the lower energy state due to greater overlap with the gain region from the Gaussian pumps. In contrast to Fig. S10a the flatband between the X and the M points in Fig. S10b has now obtained the highest gain due to this non-Hermitian splitting between the x, y P -orbitals in the A and C sublattices.