Nonlinear Polariton Fluids in a Flatband Reveal Discrete Gap Solitons

TL;DR

This paper experimentally investigates nonlinear polariton dynamics in a flat band Lieb lattice, revealing discrete gap solitons with unique properties and complex multistability, advancing understanding of flat band physics and potential photon phase control.

Contribution

It demonstrates the formation of discrete gap solitons in a driven-dissipative flat band system, a phenomenon previously observed only in closed systems, highlighting new nonlinear behaviors.

Findings

Observation of quantized-size gap solitons

Detection of abrupt edges indicating frozen propagation

Identification of multistability in nonlinear domains

Abstract

Phase frustration in periodic lattices is responsible for the formation of dispersionless flat bands. The absence of any kinetic energy scale makes flat band physics critically sensitive to perturbations and interactions. We report here on the experimental investigation of the nonlinear dynamics of cavity polaritons in the gapped flat band of a one-dimensional Lieb lattice. We observe the formation of gap solitons with quantized size and very abrupt edges, signature of the frozen propagation of switching fronts. This type of gap solitons belongs to the class of truncated Bloch waves, and had only been observed in closed systems up to now. Here the driven-dissipative character of the system gives rise to a complex multistability of the nonlinear domains generated in the flat band. These results open up interesting perspective regarding more complex 2D lattices and the generation of correlated photon phases.

Figures

Equations

$$\displaystyle\hat{H}=\sum_{l,n}E_{l}|l_{n}\rangle\langle l_{n}|-\sum_{n}\Big{(}t\big{(}|B_{n}\rangle\langle C_{n}|$$

$$\displaystyle+|B_{n}\rangle\langle C_{n+1}|\big{)}+t^{\prime}|A_{n}\rangle\langle B_{n}|+\mathrm{h.c.}\Big{)}$$

$$i\hslash \frac{\mathrm{d}{\psi }_n(t)}{\mathrm{d}t}=$$

$$-\sum _m{t}_{nm}{\psi }_m(t)+iF{e}^{-x_n^2/4{\sigma }^2}{e}^{-i{k}_p{x}_n}$$

$$t_{nm}=\{\begin{array}{ll}t& \mathrm{if}\mathrm{sites}n,m\mathrm{are}\mathrm{neighboring}B\mathrm{and}C\mathrm{sites},\\ t^{\prime }& \mathrm{if}\mathrm{sites}n,m\mathrm{are}\mathrm{neighboring}A\mathrm{and}B\mathrm{sites},\\ 0& \mathrm{otherwise},i.e.\mathrm{if}\mathrm{sites}n,m\mathrm{are}\mathrm{not}neighbors.\end{array}$$

Full text

Geometric frustration in periodic quantum media is responsible for the existence of energy bands with no dispersion. These systems extremely sensitive to any perturbation like disorder or interaction. Fascinating many body physics occurs in flatbands, among which the formation of spin liquids and spin ices Balents (2010), itinerant ferromagnetism Lieb (1989); Tasaki (1992), fractional quantum Hall states Tsui et al. (1982) or superconductivity in twisted bilayer graphene Cao et al. (2018a, b). In the quest for emulation of this rich phenomenology in controlled systems, pioneering works on frustrated lattices with a flatband have been realized in recent years Leykam et al. (2018), using different analog systems like cold atoms Jo et al. (2012); Taie et al. (2015), arrays of coupled optical waveguides Guzmán-Silva et al. (2014); Mukherjee et al. (2015); Vicencio et al. (2015) or semiconductor microcavities Masumoto et al. (2012); Jacqmin et al. (2014); Baboux et al. (2016); Klembt et al. (2017); Whittaker et al. (2018). However, most of these works have been limited so far to the linear regime, where particle-particle interactions are negligible. Despite interesting theoretical predictions Vicencio and Johansson (2013); Di Liberto et al. (2018), the experimental exploration of many-body physics in synthetic frustrated flatbands remains at its infancy.

Exciton-polaritons in semiconductor microcavities have emerged as a powerful platform to study quantum fluids in a driven-dissipative context Carusotto and Ciuti (2013). Polaritons are quasi-particles arising from the strong radiative coupling between photons confined in a semiconductor microcavity and excitons confined in quantum wells (QW). Their excitonic component provides significant repulsive interactions (equivalent to a Kerr-like nonlinearity), and the dissipative nature of the system allows for direct injection of polaritons at a given energy. This has enabled the observation of hydrodynamic features including superfluidity Amo et al. (2009), nucleation of vortices Nardin et al. (2011); Sanvitto et al. (2011) and solitons Grosso et al. (2011); Amo et al. (2011). Additionally, the potential landscape probed by polaritons can be sculpted using a variety of techniques Schneider et al. (2017). This allows controlling the polariton band structure by engineering lattices and offers a versatile playground to investigate the interplay between kinetic and interaction energy. Recently several groups have reported the realization of a flatband for polaritons: fragmentation of a polariton condensates induced by disorder was observed Baboux et al. (2016), as well as interesting polarization textures induced by spin-orbit coupling Whittaker et al. (2018).

In this Letter, we report the experimental investigation of the nonlinear response of polaritons in the flatband of a one-dimensional Lieb lattice of coupled micropillar cavities. We observe the generation of bright quantized nonlinear domains with abrupt and well-defined edges. These sharp profiles are due to the fact that the interaction energy induced by the drive cannot be accommodated as kinetic energy in the flatband. As a result, propagation of switching fronts is frozen. The size of the domains evolves through abrupt jumps as the pumping power is swept, and multistability is evidenced around each of these jumps. Theoretical analysis of the observed features indicates that these domains belong to the family of gap solitons named truncated Bloch waves Alexander et al. (2006); Wang et al. (2009); Anker et al. (2005); Bennet et al. (2011); Bersch et al. (2012), which had never been observed in a driven-dissipative context.

The 1D Lieb lattice (see Fig. 1(a)) is one of the simplest lattices hosting a flat energy band Lieb (1989). The unit cell (UC) contains three -$A,B,C$- sites linked by nearest-neighbor couplings $t$ ($t^{\prime}$) between $B$ and $C$ ($A$ and $B$) sites. In the tight-binding approximation, the corresponding single-particle Hamiltonian is:

[TABLE]

where $|l_{n}\rangle$, with $l\in\{A,B,C\}$, is the state on site $l$ in the $n$th unit cell, with on-site energy $E_{l}$. The energy spectrum of the 1D Lieb lattice presents three bands. While eigenfunctions in the lower (resp. upper) band are constructed with the same (opposite) phase for neighboring sites, eigenfunctions for the middle band show alternating phase sign on the $A,C$ sublattice. This creates a phase frustration on sites $B$: when $E_{A}=E_{C}$, a destructive interference induces zero wave function amplitude on $B$ sites. As a result, the middle band is flat with localized eigenstates of the form $|f_{n}\rangle=\left(|A_{n-1}\rangle-|C_{n}\rangle+|A_{n}\rangle\right)/\sqrt{3}$, as shown in Fig. 1(a). When $E_{A}\neq E_{C}$, the interference on $B$ sites is not fully destructive and the middle band becomes dispersive.

To experimentally implement a polariton Lieb lattice, we use a semiconductor heterostructure grown by molecular beam epitaxy. It consists of a $\lambda$ GaAs layer embedded between two $\mathrm{Ga_{0.9}Al_{0.1}As/Ga_{0.05}Al_{0.95}As}$ distributed Bragg reflectors (DBR) with 36 (top) and 40 (bottom) pairs (quality factor $Q\approx 50000$). A single 8 nm $\mathrm{In_{0.05}Ga_{0.95}As}$ QW is inserted at the center of the cavity, resulting in a $3.5\ \mathrm{meV}$ Rabi splitting. The cavity is processed into an array of coupled micropillars (Fig. 1(b)) using electron beam lithography and dry etching down to the GaAs substrate. Each micropillar is mapped to a site of the tight-binding Hamiltonian while the coupling between sites is provided by the finite overlap between neighboring pillars. The on-site energy of the lowest energy mode in each micropillar is controlled by the pillar diameter, but also depends on the number of adjacent pillars, whose proximity reduces the local confinement. Thus a fine tuning in the pillar diameters is required to match $A$ and $C$ on-site energies and obtain a flatband. We choose a 2.4 $\mathrm{\mu m}$ distance between adjacent pillars (resulting in UC size $a=4.8\ \mathrm{\mu m}$), and $3.0\,\mathrm{\mu m}$, $2.8\,\mathrm{\mu m}$ and $2.9\,\mathrm{\mu m}$ for the diameters of $A$, $B$ and $C$ pillars.

Optical spectroscopy is performed at 4 K, using a tunable continuous wave (cw) monomode excitation laser focused on the 1D lattice. The signal is collected in transmission geometry through the back of the sample, using a lens with 0.5 numerical aperture, and focused on the entrance slit of a spectrometer coupled to a CCD camera. The emission is analyzed in real- or momentum-space, by imaging either the sample surface or the Fourier plane of the collection lens. Using a $\lambda$/2 wave-plate and a polarizer, we select the linearly polarized emission either parallel (H), or orthogonal (V) to the lattice.

First, we characterize the band structure of the polariton lattice in the linear regime. The lattice is excited non-resonantly tuning the laser energy around 1.6 eV, with a Gaussian-shaped elongated spot and weak pumping power. The energy difference between the uncoupled cavity mode and the exciton resonance amounts to $-3.5\ \mathrm{meV}$. The momentum-space resolved emission for both H and V polarization is shown in Fig. 1(c). Three bands are evidenced arising from hybridization of micropillar confined modes. They are well reproduced by tight-binding calculations, except for the upper band where the expected band folding is not observed in the experiment. This deviation can be explained in terms of mixing with higher-energy bands, neglected by the tight-binding description. For H polarization, the middle band is dispersionless and is gapped from the two other bands. Because of polarization dependent boundary conditions for the electromagnetic field, the flatband condition cannot be achieved simultaneously for two orthogonal linear polarizations Baboux et al. (2016). Indeed, under V polarization all bands are dispersive. The faint intensity modulation visible in Fig. 1(c) arises from multiple reflections between the bottom mirror and the polished back side of the substrate.

To investigate the polariton nonlinear response in the flatband, we inject polaritons with a quasi-resonant laser at energy detuning $\Delta$ from the flatband. In the experiments described in the following, the flatband is at energy $E_{C}=1468.9\ \mathrm{meV}$ and the laser energy $\hbar\omega_{H}$, is blueshifted by $\Delta=90~{}\mathrm{\mu eV}$, corresponding to about one third of the energy gap separating the flatband from the upper one (Fig. 1(c)). The pumping beam is polarized along H and focused into a Gaussian-shaped elongated spot, of $40\,\mathrm{\mu m}\times 3\,\mathrm{\mu m}$ FWHM (Fig. 2(c)). Its $5.0^{\circ}$ angle of incidence matches the edge of the first Brillouin zone (BZ), i.e. $k=\pi/a$, enabling efficient coupling to the flatband modes Sup .

Figure 2(a) presents the measured total emitted intensity when increasing the excitation power $P$. Several abrupt jumps are observed separated by plateaus where the intensity weakly varies. To discuss the origin of the observed features, we show in Fig. 2(d-g) spatial emission patterns measured for different excitation powers. For $P=7~{}\mathrm{mW}$, above the lowest power intensity jump, Fig. 2(d) evidences the formation of a 4 UCs nonlinear domain. It is located around the center of the excitation spot and its shape does not evolve when $P$ is further increased, up to a power of 10 mW, where the next jump happens. We then observe the formation of a larger nonlinear domain with 10 UCs. Actually, every jump corresponds to an evolution of the domain size by a discrete number of UCs, as summarized in Fig. 2(h). The size of the domains is well defined because their edges are extremely sharp: the emission intensity drops by more than an order of magnitude over one UC on each side.

To get more insight into the physics, we solve the steady-state of a discretized Gross-Pitaevskii equation that includes pump and loss terms Carusotto and Ciuti (2004). The evolution of the polariton amplitude $\psi_{n}$ on site $n$, under a cw resonant drive, is governed, in the frame rotating at the drive frequency $\omega$, by the equation:

[TABLE]

where $\hbar g$ is the polariton-polariton interaction constant, $\gamma$ the polariton linewidth. $E_{n}$ is the on-site energy and $t_{nm}$ are the couplings to neighboring sites, deduced from $t$ and $t^{\prime}$ defined before Sup . $F$ is the drive amplitude, $k_{p}=\pi/a$ the mean wavevector of the drive and $x_{n}$ the spatial position of site $n$. The spatial distribution of the drive excitation is chosen Gaussian with $\sigma=3.5$ UCs and its detuning from the flatband, defined as $\Delta=\hbar\omega-E_{C}$, is $\Delta=3\gamma=90\ \mathrm{\mu eV}$.

Figure 2(b) presents the calculated total intensity in a 40 UC chain when increasing the drive amplitude. A series of abrupt jumps is observed in good qualitative agreement with the experiment. Similar to the experiment, the simulated steady-state spatial profiles consist in nonlinear domains of finite size Sup and each jump in the total intensity corresponds to a discrete change in the domain size.

The origin of these quantized nonlinear domains in the flatband is intimately linked to the fact that the pump energy lies in an energy gap. Locally, the system switches into the nonlinear regime (with high polariton intensity) when the local pump intensity is higher than a threshold value so that the local interaction energy $\hbar g|\psi|^{2}\gtrsim\Delta$. Outside this high excitation region, the system remains in the linear regime with very low intensity because the excitation laser lies within an energy gap. This is the reason why abrupt edges are formed. The size of the domain is determined by the number of UCs for which the local pump intensity is higher than the threshold value. Because of the Gaussian shape of the excitation spot, the threshold condition is reached for more and more UCs as the power is increased, and domains of increasing size are formed. If the excitation profile was completely flat, all excited UCs would switch simultaneously.

These nonlinear states belong to the general family of gap solitons. They have been discussed in a different context Alexander et al. (2006); Wang et al. (2009) and named truncated Bloch waves (TBW), because their pattern is similar to a spatial portion of the excited Bloch states. Notice that in our system, TBW can also be observed when exciting in a gap above a dispersive band, as illustrated by simulations presented in the supplemental Sup . In a flatband, there is no kinetic energy to overcome so that discrete TBW are formed as soon as the laser detuning overcomes the band linewidth.

A significant difference between numerical simulations and experiments is revealed in the quantitative domain size versus $F$. Once the first domain is formed (at $F^{2}=120$), the calculated domain size increases by exactly one UC at each jump (see Fig. 2(i)). The situation is different in the experiment where for instance an increase of 6 UCs is observed at 10 mW (Fig. 2(h)). In the supplemental Sup , we show that disorder in the on-site energies explains this apparent discrepancy. Experimentally, disorder mainly stems from small fluctuations in the pillar size and shape, caused by etching. A local redshift of the flatband eigenstates acts as a barrier for the nonlinear domain: when the excitation strength is sufficient to overcome this barrier, the domain size can abruptly increase by several unit cells, as observed experimentally. Notice also that the drive with wave vector at the edge of the BZ imposes a $\pi$ phase difference between adjacent UCs. This creates destructive interference on $A$ sites explaining why $A$ sites are generally dark within the nonlinear domains. Bright $A$ sites are only observed at the edges of the nonlinear domains or in regions where disorder is of the order of the interaction energy, thus altering the interference (see for instance Fig. 2(d) or Fig. 2(e), left side).

TBW have been experimentally observed in closed systems Anker et al. (2005); Bennet et al. (2011); Bersch et al. (2012), but to our knowledge have never been reported in a driven-dissipative system. Associated to the presence of dissipation, such nonlinear photonic systems are expected to present hysteretic behaviors Carusotto and Ciuti (2013); Baas et al. (2004). To probe this property, we excite the flatband with $\Delta=90\ \mathrm{\mu eV}$, and scan the power up and down around each nonlinear jump of the intensity. This reveals multiple hysteresis cycles, summarized in Fig. 3(a). In fact, each abrupt jump is associated with a hysteresis cycle. Depending on the value of the power and on the power ramp history, several configurations of nonlinear domains can be achieved. As an illustration, Fig. 3(b-e) show four nonlinear spatial patterns which can be generated for $P=10.2$ mW. Each of them shows a TBW soliton, with well defined number of bright unit cells and no emission outside the nonlinear domains. Thus the polariton nonlinear response of the flatband reveals complex multistable behavior when scanning the power up and down.

We now compare the experiments described above with the nonlinear response when pumping the system within a dispersive band. Direct comparison can be done in the very same lattice simply by injecting polaritons with an excitation linearly polarized along the V direction. The polarization dependent on-site energies result in a middle band that is dispersive in this case (Fig. 1(c), right panel). We use pumping conditions comparable to those used in Fig. 2: $\Delta=60~{}\mathrm{\mu eV}$ with respect to the bottom of the band, and same angle of incidence (Fig. 1(c)). The results are presented in Fig. 4(a): a nonlinear increase in the total emitted intensity is observed at $P=10~{}\mathrm{mW}$, followed by an abrupt intensity jump at $P=17~{}\mathrm{mW}$, associated with a hysteresis cycle. The emission spatial pattern measured in the nonlinear regime for $P=19~{}\mathrm{mW}$ is shown in logarithmic scale in Fig. 4(c). Polariton propagation out of the pumping region is clearly evidenced: the emission spreads over the entire portion of the lattice under investigation, in stark contrast with the case of the flatband (also shown in logarithmic scale in Fig. 4(d)). Exponential fits to the intensity profiles (see Fig. 4(e,f)) allow estimating a polariton propagation distance of $13.6~{}\mathrm{\mu m}$ (2.8 UCs) for the dispersive band, to be compared to only $2.1~{}\mathrm{\mu m}$ (0.4 UCs) for the flatband. Polariton propagation outside the high density region is possible in the dispersive band thanks to the conversion of interaction energy into kinetic energy. The resulting propagation of switching fronts has been previously used to realize a spin switch Amo et al. (2010). When injecting the quantum fluid into a gap, such propagation is totally suppressed.

We simulated the resonant excitation experiment for the case of the dispersive band. The calculated total intensity shown in Fig. 4(b) presents a single hysteresis cycle, in agreement with previous reports Carusotto and Ciuti (2013); Baas et al. (2004). In the supplemental, we show that the first nonlinear increase observed in the experiment at $P=10~{}\mathrm{mW}$ is well accounted for by introducing disorder in the simulation Sup . The spatial profiles calculated for the dispersive and the flatband with parameters corresponding to Fig. 4(e,f) (no disorder) are shown in Fig. 4(g,h). They nicely reproduce the difference in intensity profile when propagation of the switching front is frozen (flatband) or not (dispersive band).

In conclusion, we have shown that the nonlinear response of a polariton fluid resonantly injected into a flatband is governed by the emergence of gap solitons of the family of TBW. They are discrete nonlinear domains whose abrupt edges reflect the freezing of kinetic energy, and show complex multistable patterns under non-homogeneous spatial excitations. Similar experiments in 2D lattices with geometric frustration, where the flatband is touching a dispersive band (kagome) or crosses Dirac points (2D Lieb lattice), could offer interesting perspectives, as well as the investigation of the Bogoliubov excitations when the flatband is driven in the nonlinear regime. Using polariton structures with stronger interaction strength Togan et al. (2018); Rosenberg et al. (2018); Tsintzos et al. (2018); Knüppel et al. (2019), experiments beyond semi-classical approximation could be envisioned: in this context, flatbands are also particularly relevant for the study of many-body correlated phases Wu et al. (2007); Sun et al. (2011); Neupert et al. (2011); Clark et al. (2019).

.1 Supplemental Material: Nonlinear Polariton Fluids in a Flatband Reveal Discrete Gap Solitons

I Resonant excitation of the flatband modes

To determine the optimal excitation scheme for the resonant drive of the flatband, we study the momentum-space resolved photoluminescence (PL) of the flatband modes. This pattern can be obtained by spectrally filtering the emission, measured under weak non-resonant pumping, at the flatband energy. The $(k_{x},k_{y})$ map of the emission is then reconstructed from spectra measured at different values of $k_{y}$, such as the one shown in Fig. 1(c) of the main text. The result is presented in Fig. S1(b): the intensity is zero at the center of the Brillouin zone (BZ) ($(k_{x},k_{y})=0$), and is maximal at the BZ edges ($k_{x}=\pi/a$). This reflects the antisymmetric nature of the flatband eigenmodes (opposite phase on $A,C$ pillars). Thus, to ensure efficient coupling to the polariton states in the flatband, the pumping beam is tilted from normal incidence by $5.0^{\circ}$ along the $x$ direction, corresponding to the edge of the first BZ.

The real-space PL at the energy of the flatband can also be obtained with a similar method (spectral filtering of the real-space PL) and is shown at the top of Fig. S1(b). As expected due to geometric frustration, we measure vanishing intensity on $B$ sites.

Let us comment briefly on the consequence of exciting the flatband at the BZ edge. When the wave-vector of the driving field is equal to $k=\pi/a$, it imposes a phase difference of exactly $\pi$ between neighboring unit cells. The opposite sign on $C$ sites in neighboring unit cells leads to a destructive interference on $A$ sites. As an illustration, Fig. S1(c,d) shows calculated real- and momentum-space emission pattern of two different localized eigenstates: a single plaquette $|f_{n}\rangle$, and a linear superposition of four plaquettes, of same magnitude but alternating sign on neighboring unit cells, which can be written $\sum_{j=0}^{3}(-1)^{j}|f_{n+j}\rangle$. To compute these radiation patterns, we use a simplistic description of the eigenfunctions of the chain of pillars: we consider a Gaussian-shaped orbital per pillar (corresponding to the $s$ mode). To construct a given wave function, we assign the amplitude and phase computed from the tight-binding model to each of these Gaussian shaped orbitals. The momentum-space radiation pattern is obtained by Fourier transformation of this wave function. The state in Fig. S1(d), with a phase difference of $\pi$ between neighboring plaquettes, corresponds to a Bloch state $|\psi(\pi/a)\rangle$, but truncated to only 4 unit cells.

Importantly the non-linear domains measured in our experiment indeed present a spatial pattern similar to the one calculated in Fig. S1(d). This pattern reflects the phase imposed by the drive at the edge of the BZ, resulting in low intensity on $A$ sites inside the non-linear domains because of destructive interferences. $A$ sites have significant intensity only at the edge of the domains or in regions where disorder overcomes interaction energy (like in Fig. 2(d) or Fig. 2(e) in the main text) and locally breaks the destructive interferences.

II Parameters for numerical simulations

The discrete Gross-Pitaevskii equation introduced to model our experiments (Eq. (2) of the main text) is a set of $3N$ equations, that describe the time evolution of the polariton amplitude on each site. $N$ is the number of unit cells in the lattice. The coupling terms $t_{n,m}$ between the different sites are linked to the coupling constant $t$ and $t^{\prime}$ of the Lieb tight-binding Hamiltonian (Eq. (1) of the main text) as follows:

[TABLE]

The values of the parameters are deduced from the tight-binding fit to the measured polariton dispersion as shown in Fig. 1(c) of the main text. In the case of the flatband, we take $E_{A}=E_{C}=E_{0}$, $E_{B}=E_{C}-10\gamma$ and $t=t^{\prime}=10\gamma$, with $\gamma=30\mathrm{\mu eV}$ (and the energy offset $E_{0}=0$ for simplicity). For the dispersive band, we use $E_{A}=E_{C}+6\gamma$, all other parameters being unchanged.

III Numerical simulations: steady-state spatial profiles

In this section we present the steady-state spatial intensity profiles calculated for the quasi-resonant injection of polaritons in the flatband, that correspond to the total intensity and domain size from Fig. 2(b,i) of the main text. The drive detuning is $\Delta=3\gamma$ and the drive wave vector $k_{p}=\pi/a$. The calculated total intensity versus drive power $F^{2}$ from Fig. 2(b) of the main text is reproduced in Fig. S2(a). Fig. S2(b-d) shows the calculated spatial profile on pillars $C$ for different values of $F^{2}$. In each case, a nonlinear domain delimited by a sharp drop in occupation at the edges is clearly identified. As $F^{2}$ is increased above the first abrupt intensity jump, each intensity jump corresponds to an increase of the domain size by exactly 1 unit cell (UC).

As a comparison, we repeat the numerical simulation with same excitation conditions ($\Delta=3\gamma$, $k_{p}=\pi/a$), but with $E_{A}=6\gamma$, such that the middle band is now dispersive and the drive frequency lies within the band. The results are presented in Fig. S2(e-h). As detailed in the main text, a single jump is observed in the calculated total intensity as the drive power $F^{2}$ is increased. Moreover, the spatial density profile in the nonlinear regime is smooth and does not evolve significantly as $F^{2}$ is increased.

IV Influence of disorder in a flatband

As explained in the main text, we experimentally observe that the size of the nonlinear domain formed in the flatband can jump by several unit cells at a time when the pumping power is increased (see Fig. 2(h) of the main text, reproduced in Fig. S3(a)). This feature is not reproduced by the simulation (see Fig. 2(i) of the main text), where increases by one unit cell only are obtained (except for the first jump). In the following we show that disorder in the on-site energies can explain this discrepancy.

Disorder can strongly affect the physics of particles in a flatband, for example leading to the fragmentation of a bosonic condensate into plaquette-sized localized modes Baboux et al. (2016). Indeed, since kinetic energy is zero in a flatband, any finite amount of disorder will break the flatband picture. In a dissipative context, disorder strength needs to be greater than the linewidth to significantly alter the physics. Experimentally, disorder mainly stems from small fluctuations in the pillar size and shape, caused by etching.

An estimate of the disorder strength can be extracted from resonant spectroscopy of the flatband eigenstates in the linear regime. Figure S3(b) shows the light intensity transmitted through $A$ and $C$ pillars when scanning the laser energy. When the laser is in resonance with an eigenstate, an intensity maximum is observed. The figure clearly indicates some spatial energy spreading of the eigenstates across the lattice. More precisely for this particular part of the chain, a redshift is observed to the left of UC 0, and to the right of UC 5. The maximal energy difference between the different states is around $80\,\mathrm{\mu eV}$, comparable to the laser detuning $\Delta=90\mathrm{\mu eV}$ used in Fig. S3(a) (and Fig. 2 of the main text). Thus in the experiments disorder strength is comparable to the interaction energy. Note that imaging of the eigenstates as done in Fig. S3(b) does not allow extracting precise on-site disorder on each individual pillar. Nevertheless, to get a flavor of the effect of disorder on the nonlinear domains, we introduce in our simulations a distribution of on-site energies which results in a distribution of eigenstates resembling the measured one. A redshift $\delta_{\mathrm{dis}}$ is introduced for the on-site energy of all sites on 2 UCs to the left of the excitation spot (see Fig. S3(c)). The corresponding simulation of the nonlinear domain size versus $F^{2}$ is shown in Fig. S3(d) for different disorder strengths. As in the experiment, we observe series of jumps of different amplitudes. For instance for $\delta_{\mathrm{dis}}=0.5\gamma$, an abrupt jump from 6 to 9 UCs occurs at $F^{2}\approx 250$. It corresponds to a progression of the domain edge through all redshifted sites at once. For stronger disorder amplitude, additional big jumps in the domain size are observed at higher excitation powers. The redshifted sites thus create a barrier for the domain edges, and modify the growth of the domains with power. In the experiment, the disorder landscape is certainly more complex but this simple simulation provides good understanding of the effect of disorder on the observed nonlinear dynamics. We have verified on several lattices realizing different disorder configurations that the nonlinear behavior reported here is qualitatively the same.

V Truncated Bloch Waves in the gap above a dispersive band

We investigate with numerical simulations the behavior of a nonlinear fluid injected in the gap above a dispersive band. Figure S4(b-d) presents the steady-state profiles calculated in the nonlinear regime and without disorder, for $E_{A}=6\gamma$ and different values of the drive energy detuning with respect to the bottom of the middle band: $\Delta=4$, $5$ and $7\gamma$. Note that for $E_{A}=6\gamma$, the width of the middle band is $\sim 4.6\gamma$, so that when $\Delta>4.6\gamma$ the drive lies within the gap. For $\Delta=4\gamma$, i.e. for a drive below the band edge, the propagation outside the spot is visible as a spatial exponential decay of the intensity. The propagation length $L$ characteristic of the spatial decay is given by $L=v_{g}/\gamma$, with $v_{g}=\hbar^{-1}(\partial E/\partial k)$ the group velocity at energy $\Delta$. Increasing the drive energy to $\Delta=5\gamma$, a sharp spatial decrease in the intensity is now observed at UC $\pm 11$. For $\Delta=7\gamma$, further into the gap, the domain edge is even sharper. In this excitation configuration, since the drive injects polaritons within the gap, there is no single-particle state at this energy. As a result, the interaction energy provided by the drive cannot be converted into kinetic energy: propagation of particles out of the excitation region is prevented. This localization mechanism arising from the interplay between interactions and the existence of an energy gap is precisely the one at play in the formation of gap solitons, and in particular of Truncated Bloch Waves, as originally discussed in Refs. Alexander et al. (2006); Alexander and Kivshar (2006); Wang et al. (2009).

Thus when the dispersive band is excited within the gap at high energy, Truncated Bloch waves are excited in a similar way as for the flatband. The energy injected in the system is larger than the maximum kinetic energy the system can accommodate so that non linear domains with sharp edges are formed. In the flatband, since kinetic energy is strictly zero, this regime is achieved as soon as the driving energy overcomes the other energy scales of the system, namely the spectral linewidth and disorder.

VI Multistability of the nonlinear domains

In Fig. S5 we present several experimental power scans obtained with excitation parameters similar as those used in Fig. 2 of the main text ($\Delta=90\ \mathrm{\mu eV}$ and $k_{p}=\pi/a$). For each of these power scans, the starting condition is denoted by a black square. Fig. 3(a) of the main text reproduces all these measurements on top of each other.

VII Influence of disorder within a dispersive band

Disorder also has an influence on the nonlinear regime in the dispersive band. This is due to the fact that the disorder amplitude in the experiment, on the order of $80\ \mathrm{\mu eV}$, is comparable to the interaction and kinetic energy of the fluid with our choice of laser detuning $\Delta=60\ \mathrm{\mu eV}$. In Fig. S6, we present the results of a numerical simulation taking into account disorder in the dispersive band: we introduce a redshift $\delta_{\mathrm{dis}}=2\gamma$ on the same sites as in Fig. S3(c), $\Delta=1.5\gamma$ and $E_{A}=6\gamma$. The total population versus $F^{2}$ in the up and down scans are in excellent agreement with the experimental results from Fig. 4(a) of the main text, reproduced in Fig. S6(b). Indeed, the presence of disorder explains the first nonlinear increase in the total intensity before the abrupt jump (only one jump was observed in disorder free simulations, see Fig. 4(b) of the main text).