# Quantum fluids of light in all-optical scatterer lattices

S. Alyatkin  Skolkovo Institute of Science and Technology, Moscow, Bolshoy
Boulevard 30, bld. 1, 121205, Russia  H. Sigurdsson  Skolkovo Institute of
Science and Technology, Moscow, Bolshoy Boulevard 30, bld. 1, 121205, Russia
A. Askitopoulos  Skolkovo Institute of Science and Technology, Moscow, Bolshoy
Boulevard 30, bld. 1, 121205, Russia  J. D. Töpfer  Skolkovo Institute of
Science and Technology, Moscow, Bolshoy Boulevard 30, bld. 1, 121205, Russia
P. G. Lagoudakis<EMAIL_ADDRESS>Skolkovo Institute of Science and
Technology, Moscow, Bolshoy Boulevard 30, bld. 1, 121205, Russia

(August 27, 2024)

###### Abstract

One of the recently established paradigms in condensed matter physics is
examining a system’s behaviour in artificially constructed potentials, giving
insight into physical phenomena of quantum fluids in hard-to-reach settings. A
prominent example is the matter-wave scatterer lattice, also known as the
barrier lattice or repulsive Dirac comb. There, high energy matter waves
undergo transmission and reflection through narrow width barriers leading to
stringent phase matching conditions with subsequent lattice band formation. It
is one of the most well taught system in quantum mechanics but its realisation
for macroscopic matter-wave fluids has remained elusive, in contrast to
evanescently coupled lattice sites or waveguides. Here, we implement and study
a system of exciton-polariton condensates in a non-Hermitian Lieb lattice of
scatterer potentials by optically injecting incoherent exciton clouds which
both emit, and interact with traveling polariton waves. By fine tuning the
lattice parameters, we reveal a nonequilibrium phase transition between two
distinct regimes of polariton condensation: a scatterer lattice of gain guided
polaritons condensing on the lattice potential maxima, and trapped polaritons
condensing in the lattice potential minima. The transition is characterised by
multimodal condensation due to gain competition between the two regimes.
Energy tomography on the polariton emission enables us to measure the
intricate band structure of the optically induced lattices. Our results pave
the way towards unexplored physics of non-Hermitian fluids in non-stationary
mixtures of confined and freely expanding waves.

## Introduction

Artificial lattices are highly attractive for obtaining insight into
properties of crystal structures in the solid state, and for creating
patterned structures not found in nature. They can be used to investigate spin
frustration in triangular geometries Skjaervo_NatRev2020 , massless Dirac
fermions in honeycomb structures Wang_NatNano2018 , and strongly correlated
states in dispersionless flat bands Polini_NatureNano2013 . Conventionally,
the physical properties of lattices are investigated through an appropriate
choice of confined single-particle states (Wannier functions), such as
electrons bound to their atoms, ultracold atoms in optical traps
Schafer_NatRev2020 , or index-guided electromagnetic waves in photonic
crystals Joannopoulos_photonic_crystals . The evanescent transfer of energy
between lattices sites is then captured within the tight-binding approximation
to explain crystal band formation, the centerpiece of solid-state physics
ashcroft_solid .

However, in contrast to lattices of tightly confined states there exists the
inverse case of coherent matter-wave scattering in the diffractive regime of
lattices made up of repulsive potentials much smaller in size than the lattice
spacing. There, waves with wavelength smaller than the lattice constant
experience strong reflection and diffraction from the lattice (e.g., Bragg’s
law). Since the early work of Kronig and Penney Kronig_PRSL1931 , an ordered
arrangement of spherically symmetric scatterers is known to give rise to bands
and bandgaps. The study of elastic scattering of incoming and outgoing states
on a static object is at the heart of mechanical, electromagnetic-,
Schrödinger-, and sound-wave scattering theory Newton_ScatteringTheory .
However, the realisation of lattices of scatterers acting on matter waves with
coherence length exceeding the scatterer spacing is much harder than the
implementation of lattices of confined (evanescently coupled) wavefunctions,
and the connection between the two regimes remains mostly unexplored.

Figure 1: Schematic of the optical excitation pattern and resulting polariton
condensates. Excitation intensity profile composed of Gaussian pump spots
arranged in (a,b) a Lieb pattern for two different lattice constants $D$ and
(c) an inverse Lieb pattern with the potential minima (dark areas) forming a
conventional Lieb lattice. (d-f) Corresponding black-white surface plots of
the pump induced potential landscapes with the polariton condensate density
overlaid as red envelope. In (d) polariton condensation occurs on the
potential maxima (gain guided condensates) resulting in highly energetic
(ballistic) condensate waves whereas in (e,f) condensation takes place in the
potential minima between pump spots.

Semiconductor microcavities in the strong coupling regime are especially
appealing for engineering artificial lattices as they host matter-wave modes
known as exciton-polaritons (from here on polaritons)
kavokin_microcavities_2007 . These modes possess large nonlinearities,
picosecond scale response times, and permit easy optical write-in and read-out
of information. Polaritons can undergo power-driven nonequilibrium Bose-
Einstein condensation, making them favourable candidates to study low
threshold room-temperature lasing Christopoulos_PRL2007 , optical many-body
hydrodynamic phenomena Stepanov_NatComm2019 ; Suarez_LSA2020 , topological
phases klembt_exciton-polariton_2018 ; Gianfrate_Nature2020 , and
implementation of optical information processing Zasedatelev_NatPho2019 .
Moreover, strong interparticle interactions result in repulsion of condensate
polaritons from a background of uncondensed particles (i.e., photoexcited
exciton reservoir), co-localised with the pumped area Wertz_NatPhys2010 ;
askitopoulos_polariton_2013 . This enables all-optical design of non-Hermitian
(i.e., optical gain and blueshift) potential landscapes for polaritons
tosi_sculpting_2012 ; alyatkin_PRL ; Pickup_NatComm2020 , in analogy to dipole
optical traps for cold atoms Bloch_NatPhys2005  or photorefractive crystals
Neshev_PRL2004 . Indeed, with non-Hermitian potentials the wavefunction norm
is no longer conserved and, in conjunction with stabilising condensate
nonlinearities, synchronisation can spontaneously appear amplifying the
matter-wave similarly to phase-locked laser arrays Glova_QuaElec2003 ;
Mawst_IEEE2003 .

Here, we demonstrate a non-resonantly optically imprinted repulsive
(scatterer) lattice (see Figs. 1a,d) wherein scattered high energy polariton
waves, emitted from the pump areas, result in robust interference patterns due
to their ability to dynamically adjust their phase in order to condense into
the highest gain Bloch state. Our scatterer lattice is chosen to have the
edge-centred square (Lieb) arrangement, a configuration not found usually in
nature, which offers comparison against the conventional tight-binding Lieb
lattice (see Figs. 1c,f) which we also all-optically engineer. We perform full
momentum-energy space tomography to unveil the engineered lattice band
structures and their reshaping by altering the lattice parameters. We observe
a gradual nonequilibrium phase transition from the scatterer lattice of
ballistically expanding polariton waves (Figs. 1a,d) to a tight-binding
lattice of trapped condensates (Figs. 1b,e) bridged by an unstable regime of
multimodal condensation due to gain competition. Moreover, underscoring the
flexibility of our optical approach, we provide direct observation of
dispersionless P-flatband condensation achieved by using an excitation profile
forming an “inverse” Lieb lattice (Figs. 1c,f) in the same spirit as the
vacancy lattice created in electronic systems Drost_NatPhys2017 ;
Slot_NatPhys2017  or optical lattices of cold atoms Shen_PRB2010 ;
Gross_Science2017 . We point out that the majority of our findings are not
strongly dependent on the choice of Lieb lattice arrangement and can be
extended to other types such as square, honeycomb, and triangle lattices.

## Results

A strain compensated 2$\lambda$ GaAs-based planar microcavity with embedded
three pairs of In0.08Ga0.92As quantum wells cilibrizzi_polariton_2014  and an
exciton-photon detuning of $-4\;\mathrm{meV}$ is held at $\approx
4\;\mathrm{K}$ in a closed-cycle helium cryostat. The non-resonant excitation
(single mode laser tuned at 1.5578 eV) is amplitude modulated at a frequency
of 5 kHz (duty cycle 1%) with an acousto-optical modulator to avoid sample
heating. A desired pump profile is shaped by a computer-controlled reflective
phase only spatial light modulator toepfer_coherence2020  and projected onto
the sample through a microscope objective (NA=0.42). We collect the real space
polariton photoluminescence (PL), directly proportional to the condensate
density, in transmission geometry and spectrally filter it from the excitation
laser. We denote the horizontal and vertical momentum space coordinates as
$\mathbf{k}=(k_{x},k_{y})$, corresponding to the spatial frequencies along
horizontal ($x$-axis) and vertical ($y$-axis) real space coordinates,
respectively.

As mentioned above, the flexibility in tuning the lattice properties through
optical excitation structuring makes a system of microcavity polaritons
appealing to access various lattice physics in a recyclable setting. The
height of the potential landscape is determined by the excitation intensity
and, in this study, constructed by an arrangement of Gaussian shaped pump
spots. On one hand, in the scatterer lattice where the pump spots have large
separation distances with multiples of the polariton wavelength, polaritons
condense on the maxima of the pumped potentials (see Figs. 1a,d) characterised
by gain guided bright centres and ballistically expanding envelopes
Pickup_NatComm2020 ; toepfer_2020 ; alyatkin_PRL ; toepfer_coherence2020 . On
the other hand, for closely spaced pump spots, they can condense in the minima
of the potential landscape (see Figs. 1b,e and 1c,f) becoming optically
trapped Cristofolini_PRL2013 ; ohadi_synchronization_2018 ;
pieczarka_arxiv2021  partly because of their strong interactions that help
them relax in energy. In the former case, the system shares analogies with
with antiguided laser arrays Mawst_IEEE2003 , whereas in the latter case with
two-dimensional (2D) electron transport Wang_NatNano2018 .

Figure 2: (a) Measured PL from a scatterer Lieb lattice of polariton
condensates with a lattice constant set to $D=20.3$ $\upmu$m excited at
$P=P_{\mathrm{thr}}$ and (b) at $P=1.2P_{\mathrm{thr}}$. Inset in (b) shows
the calculated Bloch state density for comparison. The inset in (a) shows
sublattices denoted with letters A, B, C forming the lattice. (c,d) Energy-
resolved momentum space polariton PL at $P=1.2P_{\mathrm{thr}}$ for $k_{x}$=0
and $k_{x}=2\pi/D$ values respectively, marked on (e) with vertical lines and
yellow letters “c” and “d”. (e) Energy integrated momentum space polariton PL
at $P=1.2P_{\mathrm{thr}}$, and (f,g,h) energy-resolved “slices” of momentum
space at energies 1.6 meV, 1.3 meV and 1 meV with respect to the bottom of the
lower polariton branch (dispersion). (i) Calculated lattice bands with the red
circle corresponding to the Bloch state in the inset in (b). (j,k,l)
Numerically calculated “slices” of polariton momentum space PL from Monte-
Carlo simulations on the 2D dissipative Schrödinger equation.

In Fig. 2a we show the real space PL of an optically pumped scatterer lattice
with 96 polariton condensates at threshold power ($P=P_{\text{thr}}$) arranged
into a Lieb geometry (corresponding to Figs. 1a,d). Each condensate is
strongly gain guided with a bright centre co-localised with its respective
pump spot. The size of the Gaussian pump spots sustaining the condensates is
$\approx 2$ $\upmu$m full width at half maximum and the lattice constant is
set to $D=20.3\ \upmu$m. As schematically shown in the inset of Fig. 2a, the
Lieb lattice is composed of three square sublattices denoted with the letters
A, B, and C. Being close to threshold, the spatial coherence of each
condensate does not extend outside its respective pump spot Roumpos_PNAS2012
making them isolated (uncoupled) objects. Figures 2b and 2e show the real and
momentum space polariton PL above threshold ($P=1.2P_{\text{thr}}$), while
Figs. 2c and 2d show energy-momentum space PL along $k_{x}=0$ and
$k_{x}=2\pi/D$, respectively. By driving the system above threshold we
increase both the coherence and particle outflow from each condensate
resulting in stronger coupling between the radiating condensates and the
appearance of interference fringes (signature of synchronisation) as well as
distinct Bragg diffraction peaks in momentum space. Specifically, we observe
the condensates populating an excited Bloch state characterised by a
constructive interference peak in the centre of each lattice cell, also
visible from the calculated Bloch state shown in the inset in Fig. 2b. By
scanning the Fourier space PL (collected in reflection geometry) image across
a slit of a spectrometer with a motorised mount, we reconstruct energy-
resolved “slices” (i.e., isoenergy planes) of momentum space PL at three
specific energies, shown in Figs. 2f-h. The applied energy tomography reveals
the full picture of complex band formation (see Supplementary Movie 1). We
also calculate the band structure using two different modeling methods:
Bloch’s theorem, and Monte-Carlo sampling of the 2D dissipative polariton
Schrödinger equation [see Secs. S2 and S3 in the Supplemental Material (SM)
for details]. The Bloch analysis, as shown in Fig. 2i, reveals a zoo of
overlapping bands of distinct shapes, yet the polaritons, being so
interactive, are still easily able to relax into a definitive Bloch state
corresponding to the optimum gain (marked with a red circle and plotted in
Fig. 2b inset), in a similar spirit to coupled laser systems. The numerical
Monte-Carlo sampling of the dominant Fourier components of the scatterer
lattice shown in Figs. 2j-l gives good agreement with the experiment. We note
that the illuminated, clearly formed, bands in Figs. 2c and 2d show that
polaritons in the repulsive scatterer Lieb lattice indeed experience crystal
scattering within their coherence time. We also observed band-structure
formation for the square scatterer lattice of 5$\times$5 polariton condensates
with a similar lattice constant (see Sec. S1 in SM).

Next, we decrease the lattice constant from $D=20$ $\upmu$m to $D=10.3$
$\upmu$m and characterise the change in the polariton system behaviour.
Corresponding images with experimentally measured real space polariton PL
above condensation threshold are depicted in Figs. 3a-e. As shown in Fig. 3a,
for $D=20$ $\upmu$m, most of the PL intensity comes from the pumped areas with
clear interference fringes in time-averaged measurements as a consequence of
the scattered polariton waves leading to robust synchronisation between the
condensates toepfer_2020 ; toepfer_coherence2020 . Multiple weakly populated
energy branches collapse into one dominant mode as the pump power exceeds
threshold value, as clearly seen in Fig. 3f. Decreasing the lattice constant
to $D=16.9$ $\upmu$m results in substantial growth of PL intensity inside each
lattice cell (i.e., where a Gaussian pump spot is absent) with simultaneous
decreased PL at the pumps positions (see Fig. 3b). Driving the system above
threshold leads to dual mode condensation as confirmed by the measured
spectrum power scan shown in Fig. 3g. Decreasing the lattice constant to
$D=15.2$ $\upmu$m again dramatically modifies the polariton PL pattern. As
shown in Fig. 3c, polaritons are repelled even stronger outside the pumped
areas leading to complex PL distribution in real space. The condensate is here
fractured into multiple energy modes above threshold (see Fig. 3h), with
similar gain. Finally, decreasing to even smaller lattice constants of
$D=12.3$ $\upmu$m and $D=10.3$ $\upmu$m results in formation of trapped
condensates ohadi_synchronization_2018  as shown in Figs. 3d and 3e,
characterised by a dominant single energy mode above threshold (see Figs. 3i
and 3j). This regime is schematically depicted in Figs. 1b and 1e.

Figure 3: Optical lattices of polariton condensates pumped with Lieb geometry
for different lattice constants $D=$ 20.0, 16.9, 15.2, 12.3, 10.3 $\upmu$m.
(a-e) show corresponding real space polariton PL above condensation threshold
and (f-j) corresponding spectra as functions of pump power (energy is scaled
with respect to the bottom of lower polariton branch). Semi-transparent insets
in (a-c,e) marked with white dashed squares show results of time-averaged
numerical simulations of the condensate dynamics using the generalised Gross-
Pitaevskii equation. Inset in (d) is calculated using Bloch’s theorem. White
circles and squares overlaid with real space PL below threshold for $D=20\
\upmu$m in (k) denote spatial integration areas on top of and outside the pump
spots used to extract the contrast between gain guided ($I_{\text{spots}}$)
and trapped ($I_{\text{trapped}}$) polariton PL. This contrast is plotted as
normalised heatmaps
$S=(I_{\text{spots}}-I_{\text{trapped}})/(I_{\text{spots}}+I_{\text{trapped}})$
in (l,m) as a function of lattice constant and pump power from experiment and
simulations, respectively, revealing a gradual transition between the two
regimes. Scale bar in (a) applies also to (b-e) and (k). Yellow circles in (l)
denote the lattice constants realized experimentally.

The most intriguing physics lies in the intermediate regime where multimodal
condensation takes place between the regimes of gain guided ballistically
coupled condensates (Fig. 3a) and trapped condensates (Figs. 3d and 3e). In
order to distinguish these regimes and identify the transition between them we
introduce a contrast parameter
$S=(I_{\text{spots}}-I_{\text{trapped}})/(I_{\text{spots}}+I_{\text{trapped}})$
describing the normalised difference between the average polariton PL at the
pumping spots $I_{\text{spots}}$ (white circles in Fig. 3k) and outside the
spots $I_{\text{trapped}}$ (white squares in Fig. 3k). Here, “average” refers
to area-integrated PL divided by the integration area. The side length of the
white squares is chosen as $L=0.54D$ (the precise value does not affect the
findings) whereas the diameter $d=3.6\;\mathrm{\upmu m}$ of the white circles
is fixed for all lattice constants.

In case of well separated pump spots (i.e., $D>17$ $\upmu$m) we observe values
$S>0$ in Fig. 3l across all pump powers corresponding to high-energy gain
guided polaritons that scatter across the lattice when emitted from their pump
spots. In the opposite case of small separation distances between the pumps
(i.e., $D<12.5$ $\upmu$m) we observe $S<0$ corresponding to polaritons
becoming trapped in the lattice potential minima with dominant PL intensity
coming from the white squares in Fig. 3k. Here, high energy scattering of
polaritons across multiple lattice cells is reduced, replaced with evanescent
transfer of energy (i.e., tunneling). Therefore, at some lattice constant
(where $S\approx 0$) the polariton system undergoes a nonequilibrium phase
transition between these different condensation regimes. During this
transition, competing lattice modes fight over the gain which leads to a
“frustrated“ condensation pattern, like shown in Figs. 3b and 3c. This
transition is unique to polariton systems since it is the excitonic part of
polaritons which allows them to interact with each other and the uncondensed
exciton reservoir, which facilitates energy relaxation Wertz_NatPhys2010 . In
contrast, purely photonic systems would generally lase in their pumped gain
region. The gradual onset of trapped polaritons as the lattice constant
reduces and pump power increases (i.e., the background uncondensed reservoir
increases) can therefore be attributed to modes between the pump spots moving
into resonance with the lattice gain bandwidth due to enhanced energy
relaxation of polaritons Wouters_PRB2010 . We point out that our all-optical
experiment and extraction of $S$ for different lattice constants and pump
powers is not possible to replicate with lithographically written photonic
periodic structures where the pumped gain region would correspond to the
lattice potential minima klembt_exciton-polariton_2018 . In our experiment,
the pump gain follows the lattice potential maxima. We qualitatively reproduce
our findings through numerical mean-field simulations using the generalised
Gross-Pitaevskii equation Wouters_PRL2007  including an energy relaxation
mechanism Wouters_PRB2010  (see Sec. S4 in SM for details), shown in the
insets of Fig. 3a-c,e and Fig. 3m. The state in the inset in Fig. 3d was
calculated using Bloch’s theorem and did not appear in stable form in Gross-
Pitaevskii simulations for our chosen (fixed) set of simulation parameters.
Some discrepancy can be observed between simulation and experiment in Figs. 3l
and 3m at low powers and small lattice constants where the simulation
overestimates the gain guided polaritons. This discrepancy could be reduced by
applying a stochastic treatment (e.g., Wiener noise) to the condensate
equations of motion which would smear out the simulated condensate PL at low
powers close to threshold, or by including exciton diffusion in the model.

Figure 4: Optical inverse Lieb lattice of Gaussians with flatband
condensation. (a,e) Spatial pumping profile of Gaussian spots arranged in the
inverse Lieb geometry with lattice constants $D=13$ $\upmu$m and $16$
$\upmu$m, respectively. In (a) the cube-centred spots are 40% weaker in
intensity compared to the rest of the pump spots. In (e) all spots are equally
intense. (b,f) Real and (c,g) momentum space polariton PL above condensation
threshold corresponding to populated S-band ground state and P-flatband state,
respectively. White circles in (b) and (f) denote the pumps positions. (d,h)
Corresponding calculated real space wavefunction densities intensities from
Bloch’s theorem. Coloured circles in (a) and (h) denote the Lieb lattice unit
cells. (i) Calculated P-bands from Bloch’s theorem where the red circle marks
the flatbands. (j) Calculated dispersion cross-section at $k_{x}=2\pi/D$ from
Monte-Carlo Schrödinger numerics. Both (i) and (j) use the pump profile shown
in (e). (k,l) Calculated and measured PL from the isoenergy plane in momentum
space corresponding to the P-flatband marked with the yellow letter k in (j).
The band marked with the yellow letter k’ is the lowest energy P-band which
only appears flat between the $\Gamma$-$X$ points in the Brillouin zone.

To further demonstrate the versatility of our all-optical approach we move
away from the scatterer Lieb lattice and design now the conventional Lieb
lattice (see Figs. 1c and 1f) which has been studied vigorously in various
systems of electrons Drost_NatPhys2017 ; Slot_NatPhys2017 , photons
Guzman_NJP2014 ; Mukherjee_PRL2015 ; Vicencio_PRL2015 , cold atoms
Shen_PRB2010 ; Taie_SciAdv2015 , plasmon polaritons Kajiwara_PRB2016 , and
exciton-polaritons Klembt_APL2017 ; Whittakes_PRL2018 ; Goblot_PRL2019 ;
Harder_PRB2020 ; Scafirimuto_CommPhys2021 . It hosts an excess of phenomena
including topological phases, dispersionless flat bands, and Dirac points
making it a popular testbed in solid-state physics. In Fig. 4 we realise the
conventional Lieb lattice by arranging the Gaussian excitation spots in an
inverse Lieb lattice (Figs. 4a and 4e). Just as in Fig. 3d and 3e, here the
polariton condensates are designed to be trapped in the lattice potential
minima, interacting via tunnelling. Figures 4b and 4c show the real- and
momentum space polariton PL above threshold, respectively, in a $D=13$
$\upmu$m lattice revealing condensation into the S-band ground state in the
$\Gamma$ point at $k=0$. By adjusting the lattice parameters, the condensate
can be forced to populate a different Bloch state. Figures 4e-g show the same
measurements but for $D=16$ $\upmu$m where we now observe condensation into
P-orbitals at sublattices A and C, whereas at sublattice B destructive
interference occurs. This state belongs to a dispersionless P-flatband
Klembt_APL2017 ; Whittakes_PRL2018  and constitutes the first evidence of
polariton condensation into optically (non-Hermitian) generated flatband
states. Figures 4d and 4h show corresponding calculated lattice states using
Bloch’s theorem. The state plotted in Fig. 4h corresponds to the band denoted
by the red circle in Fig. 4i where we show the calculated P-bands along the
reduced Brillouin zone edge. We additionally calculate the dispersion of the
lattice from Monte-Carlo sampling of the dissipative Schrödinger equation and
plot the results in Fig. 4j along the $k_{x}=2\pi/D$ direction. The results
reveal bright bands where polaritons decay slowest with the P-flatband clearly
visible (marked with yellow k) and the lowest energy P-band (marked with
yellow k’) which only appears flat along the $\Gamma$-$X$ direction.

Experimentally implemented energy tomography measurements allow us to extract
PL belonging to isoenergy planes in reciprocal space which further verifies
that the condensate is populating the P-flatband. Figures 4k and 4l show the
calculated and measured PL from the P-flatband isoenergy plane, respectively,
which is marked with the yellow colored letter k in Fig. 4j. Indeed, the
similarity between Fig. 4l and the energy-integrated PL in Fig. 4g confirms
that condensate is dominantly populating a single energy state with similar
intensity pattern obtained from calculation (see Sec. S1 in SM for more
details). It should be noted that the calculated PL (Fig. 4k) shows some finer
details in momentum space which would average out in time-integrated
measurements due to non-ideal effects such as noise, cavity disorder, and pump
fluctuations. However, the qualitative structure is the same as in experiment
with most of the PL belonging to the $M$ symmetry points. The fact that the PL
is strongest at the $M$ points can be understood from the slight curvature in
the P-flatband around these points (see Fig. 4i) which facilitates the
relaxation of polaritons into the band minima.

We now address the question on why condensation occurs into the P-flatband as
shown in Figs. 4f and 4g. As mentioned in the previous paragraph, the
P-flatband in our optical lattice is not perfectly flat and is characterised
by a small dispersion (curvature) around the $M$ point which is only $\approx
5$ $\upmu$eV different in energy from both the $\Gamma$ or the $X$ points (see
Fig. 4i). This small band curvature can trigger condensation into the band
minima and is attributed to the finite potential depth of the lattice sites
(taken here to be 2 meV in calculation) which leads to deviation from the
perfectly dispersionless bands predicted by the tight binding theorem
Shen_PRB2010 . The band curvature can be reduced by increasing the potential
depth of the inverse Lieb lattice (see Sec. S5 in SM) which can be achieved by
adjusting the system properties (e.g., exciton-photon detuning or the exciton
dipole moment by appropriate choice of the semiconductor material) such that
stronger pumping (bigger blueshifting reservoir) is required to achieve
condensation. Another feature of our inverse optical Lieb lattice is that it
creates slightly elliptical confinement at sublattices A and C which splits
the energies and linewidths of the $P_{x,y}$ orbitals. In other words, the
overlap of $P_{x,y}$ orbitals into the pump (gain) region is different. For
sublattice A (C) the $P_{x}$ $(P_{y})$ orbital overlaps more with the pump
which creates a higher gain for the P-flatband state. This interpretation is
supported through a non-Hermitian tight-binding theory (see Sec. S6 in SM).

## Discussion

The demonstrated optically arranged system of exciton-polaritons presents a
flexible platform to study fundamental proposals on non-Hermitian phenomena in
artificial lattices including nonlinear reshaping of the crystal bands when
above condensation threshold. In Fig. 2 we have demonstrated the condensation
of polaritons into excited Bloch states belonging to a lattice of scatterers
(narrow waist repulsive potentials) shaped into a Lieb pattern. The observed
crystal bands and agreement with theoretical models opens up a path to explore
more intriguing effects of scattered matter-waves such as slow polaritons
Baba_NatPho2008 , guided polaritons Mekis_PRL1996 , and solitonic modes
Redondo_NatComm2014 . In Fig. 3 we have revealed a gradual nonequilibrium
phase transition from the scatterer Lieb lattice, characterised by strongly
gain guided and energetic polariton condensates, to a square lattice of
optically trapped condensates as a function of two easily tunable parameters
in experiment, the lattice constant and pump power. The transition regime is
accompanied by multimode polariton lasing which is a unique feature of
polariton systems due to their strong interactions, that leads to unexpected
condensation patterns in real space. Such multimodal behaviour of the
condensate implies close gain competition between gain guided and trapped
polaritons. Our observations therefore highlight the intriguing duality of
polaritons as “lasers” (gain guiding) and “Bose-Einstein condensates”
(thermalisation) with interesting perspectives as strongly nonlinear objects
in non-Hermitian optical lattices. As an example, the optical malleability of
our matter-wave platform allows one to study in detail the gradual departure
from a system of quantised energies to that of smoothly connected
quasimomentum states by building the lattice up cell-by-cell.

Another possible perspective is to optically engineer the lattice to probe
exotic band properties for polariton condensation. In Fig. 4 we have
demonstrated a conventional (i.e., tightly bound waves) Lieb lattice, by
packing pump spots into its inverse shape, with subsequent condensation into
flat band states. Given the nonlinear nature of exciton-polaritons,
condensation into adjustable flat band states could open a window to
investigate strongly correlated states of matter and the effects of disorder
against interactions Baboux_PRL2016 . The implemented energy tomography
methods provide complete access to the polariton states in Fourier space and
can be used to study reshaping of the spectrum under arbitrary lattice
deformation and nonlinear effects. We point out that the current study is
performed in the scalar polariton regime, but can be easily extended to
include its spin degree of freedom by changing the pump polarisation.
Furthermore, bands induced by our optically engineered landscape can be
populated using resonant excitation permitting study of the evolution of
polariton matter waves in non-Hermitian optical lattices with any chosen
crystal momentum and frequency. We believe that our work carries significant
weight in the future design and investigation of polaritonic non-Hermitian
(gain and loss) lattice physics in e.g. topological lasers Amelio_PRX2020 ;
Comaron_PRR2020 , phase transitions in many body systems Zamora_PRX2017 ;
Fruchart_Nature2021 , non-reciprocal transport Mandal_PRL2020 , and access to
a multitude of gain-induced anomalies reported for diffractive metasurfaces
Kolkowski_IEEE2020 .

## References

  * (1) Skjærvø, S. H., Marrows, C. H., Stamps, R. L. & Heyderman, L. J.  Advances in artificial spin ice.  _Nature Reviews Physics_ 2, 13–28 (2020).  URL https://doi.org/10.1038/s42254-019-0118-3.
  * (2) Wang, S. _et al._ Observation of dirac bands in artificial graphene in small-period nanopatterned gaas quantum wells.  _Nature Nanotechnology_ 13, 29–33 (2018).  URL https://doi.org/10.1038/s41565-017-0006-x.
  * (3) Polini, M., Guinea, F., Lewenstein, M., Manoharan, H. C. & Pellegrini, V.  Artificial honeycomb lattices for electrons, atoms and photons.  _Nature Nanotechnology_ 8, 625–633 (2013).  URL https://doi.org/10.1038/nnano.2013.161.
  * (4) Schäfer, F., Fukuhara, T., Sugawa, S., Takasu, Y. & Takahashi, Y.  Tools for quantum simulation with ultracold atoms in optical lattices.  _Nature Reviews Physics_ 2, 411–425 (2020).  URL https://doi.org/10.1038/s42254-020-0195-3.
  * (5) Joannopoulos, J. D., Johnson, S. G., Winn, J. N. & Meade, R. D.  _Photonic Crystals_ (Princeton University Press, 2008).
  * (6) Ashcroft, N. W. & Mermin, N. D.  _Solid State Physics_ (Brooks Cole, 1989).
  * (7) Kronig, R. D. L., Penney, W. G. & Fowler, R. H.  Quantum mechanics of electrons in crystal lattices.  _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ 130, 499–513 (1931).  URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1931.0019.
  * (8) Newton, R. G.  _Scattering Theory of Waves and Particles_ (Springer-Verlag Berlin Heidelberg, 1982).
  * (9) Kavokin, A., Baumberg, J. J., Malpuech, G. & Laussy, F. P.  _Microcavities_ (OUP Oxford, 2007).
  * (10) Christopoulos, S. _et al._ Room-temperature polariton lasing in semiconductor microcavities.  _Phys. Rev. Lett._ 98, 126405 (2007).  URL https://link.aps.org/doi/10.1103/PhysRevLett.98.126405.
  * (11) Stepanov, P. _et al._ Dispersion relation of the collective excitations in a resonantly driven polariton fluid.  _Nature Communications_ 10, 3869 (2019).  URL https://doi.org/10.1038/s41467-019-11886-3.
  * (12) Suárez-Forero, D. G. _et al._ Quantum hydrodynamics of a single particle.  _Light: Science & Applications_ 9, 85 (2020).  URL https://www.nature.com/articles/s41377-020-0324-x#citeas.
  * (13) Klembt, S. _et al._ Exciton-polariton topological insulator.  _Nature_ 562, 552–556 (2018).  URL https://doi.org/10.1038/s41586-018-0601-5.
  * (14) Gianfrate, A. _et al._ Measurement of the quantum geometric tensor and of the anomalous hall drift.  _Nature_ 578, 381–385 (2020).  URL https://doi.org/10.1038/s41586-020-1989-2.
  * (15) Zasedatelev, A. V. _et al._ A room-temperature organic polariton transistor.  _Nature Photonics_ 13, 378–383 (2019).  URL https://doi.org/10.1038/s41566-019-0392-8.
  * (16) Wertz, E. _et al._ Spontaneous formation and optical manipulation of extended polariton condensates.  _Nature Physics_ 6, 860–864 (2010).  URL https://doi.org/10.1038/nphys1750.
  * (17) Askitopoulos, A. _et al._ Polariton condensation in an optically induced two-dimensional potential.  _Phys. Rev. B_ 88, 041308 (2013).  URL https://link.aps.org/doi/10.1103/PhysRevB.88.041308.
  * (18) Tosi, G. _et al._ Sculpting oscillators with light within a nonlinear quantum fluid.  _Nature Physics_ 8, 190–194 (2012).  URL https://doi.org/10.1038/nphys2182.
  * (19) Alyatkin, S., Töpfer, J. D., Askitopoulos, A., Sigurdsson, H. & Lagoudakis, P. G.  Optical control of couplings in polariton condensate lattices.  _Phys. Rev. Lett._ 124, 207402 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevLett.124.207402.
  * (20) Pickup, L., Sigurdsson, H., Ruostekoski, J. & Lagoudakis, P. G.  Synthetic band-structure engineering in polariton crystals with non-hermitian topological phases.  _Nature Communications_ 11, 4431 (2020).  URL https://doi.org/10.1038/s41467-020-18213-1.
  * (21) Bloch, I.  Ultracold quantum gases in optical lattices.  _Nature Physics_ 1, 23–30 (2005).  URL https://doi.org/10.1038/nphys138.
  * (22) Neshev, D. N. _et al._ Observation of discrete vortex solitons in optically induced photonic lattices.  _Phys. Rev. Lett._ 92, 123903 (2004).  URL https://link.aps.org/doi/10.1103/PhysRevLett.92.123903.
  * (23) Glova, A. F.  Phase locking of optically coupled lasers.  _Quantum Electronics_ 33, 283–306 (2003).  URL https://doi.org/10.1070%2Fqe2003v033n04abeh002415.
  * (24) Mawst, L.  ”anti” up the aperture [antiguided vcsel structures].  _IEEE Circuits and Devices Magazine_ 19, 34–41 (2003).  URL https://ieeexplore.ieee.org/document/1191436.
  * (25) Drost, R., Ojanen, T., Harju, A. & Liljeroth, P.  Topological states in engineered atomic lattices.  _Nature Physics_ 13, 668–671 (2017).  URL https://doi.org/10.1038/nphys4080.
  * (26) Slot, M. R. _et al._ Experimental realization and characterization of an electronic lieb lattice.  _Nature Physics_ 13, 672–676 (2017).  URL https://doi.org/10.1038/nphys4105.
  * (27) Shen, R., Shao, L. B., Wang, B. & Xing, D. Y.  Single dirac cone with a flat band touching on line-centered-square optical lattices.  _Phys. Rev. B_ 81, 041410 (2010).  URL https://link.aps.org/doi/10.1103/PhysRevB.81.041410.
  * (28) Gross, C. & Bloch, I.  Quantum simulations with ultracold atoms in optical lattices.  _Science_ 357, 995–1001 (2017).  URL https://science.sciencemag.org/content/357/6355/995.
  * (29) Cilibrizzi, P. _et al._ Polariton condensation in a strain-compensated planar microcavity with ingaas quantum wells.  _Applied Physics Letters_ 105, 191118 (2014).  URL https://doi.org/10.1063/1.4901814.
  * (30) Töpfer, J. D. _et al._ Engineering spatial coherence in lattices of polariton condensates.  _Optica_ 8, 106–113 (2021).  URL http://www.osapublishing.org/optica/abstract.cfm?URI=optica-8-1-106.
  * (31) Töpfer, J. D., Sigurdsson, H., Pickup, L. & Lagoudakis, P. G.  Time-delay polaritonics.  _Communications Physics_ 3, 2 (2020).  URL https://doi.org/10.1038/s42005-019-0271-0.
  * (32) Cristofolini, P. _et al._ Optical superfluid phase transitions and trapping of polariton condensates.  _Phys. Rev. Lett._ 110, 186403 (2013).  URL https://link.aps.org/doi/10.1103/PhysRevLett.110.186403.
  * (33) Ohadi, H. _et al._ Synchronization crossover of polariton condensates in weakly disordered lattices.  _Phys. Rev. B_ 97, 195109 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.97.195109.
  * (34) Pieczarka, M. _et al._ Topological phase transition in an all-optical exciton-polariton lattice.  _arXiv e-prints_ arXiv:2102.01262 (2021).
  * (35) Roumpos, G. _et al._ Power-law decay of the spatial correlation function in exciton-polariton condensates.  _Proceedings of the National Academy of Sciences_ 109, 6467–6472 (2012).  URL https://www.pnas.org/content/109/17/6467.
  * (36) Wouters, M., Liew, T. C. H. & Savona, V.  Energy relaxation in one-dimensional polariton condensates.  _Phys. Rev. B_ 82, 245315 (2010).  URL https://link.aps.org/doi/10.1103/PhysRevB.82.245315.
  * (37) Wouters, M. & Carusotto, I.  Excitations in a nonequilibrium bose-einstein condensate of exciton polaritons.  _Phys. Rev. Lett._ 99, 140402 (2007).  URL https://link.aps.org/doi/10.1103/PhysRevLett.99.140402.
  * (38) Guzmán-Silva, D. _et al._ Experimental observation of bulk and edge transport in photonic lieb lattices.  _New Journal of Physics_ 16, 063061 (2014).  URL https://doi.org/10.1088%2F1367-2630%2F16%2F6%2F063061.
  * (39) Mukherjee, S. _et al._ Observation of a localized flat-band state in a photonic lieb lattice.  _Phys. Rev. Lett._ 114, 245504 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevLett.114.245504.
  * (40) Vicencio, R. A. _et al._ Observation of localized states in lieb photonic lattices.  _Phys. Rev. Lett._ 114, 245503 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevLett.114.245503.
  * (41) Taie, S. _et al._ Coherent driving and freezing of bosonic matter wave in an optical lieb lattice.  _Science Advances_ 1 (2015).  URL https://advances.sciencemag.org/content/1/10/e1500854.
  * (42) Kajiwara, S., Urade, Y., Nakata, Y., Nakanishi, T. & Kitano, M.  Observation of a nonradiative flat band for spoof surface plasmons in a metallic lieb lattice.  _Phys. Rev. B_ 93, 075126 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevB.93.075126.
  * (43) Klembt, S. _et al._ Polariton condensation in s- and p-flatbands in a two-dimensional lieb lattice.  _Applied Physics Letters_ 111, 231102 (2017).  URL https://doi.org/10.1063/1.4995385.
  * (44) Whittaker, C. E. _et al._ Exciton polaritons in a two-dimensional lieb lattice with spin-orbit coupling.  _Phys. Rev. Lett._ 120, 097401 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.097401.
  * (45) Goblot, V. _et al._ Nonlinear polariton fluids in a flatband reveal discrete gap solitons.  _Phys. Rev. Lett._ 123, 113901 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevLett.123.113901.
  * (46) Harder, T. H. _et al._ Exciton-polaritons in flatland: Controlling flatband properties in a lieb lattice.  _Phys. Rev. B_ 102, 121302 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevB.102.121302.
  * (47) Scafirimuto, F. _et al._ Tunable exciton–polariton condensation in a two-dimensional lieb lattice at room temperature.  _Communications Physics_ 4, 39 (2021).  URL https://doi.org/10.1038/s42005-021-00548-w.
  * (48) Baba, T.  Slow light in photonic crystals.  _Nature Photonics_ 2, 465–473 (2008).  URL https://doi.org/10.1038/nphoton.2008.146.
  * (49) Mekis, A. _et al._ High transmission through sharp bends in photonic crystal waveguides.  _Phys. Rev. Lett._ 77, 3787–3790 (1996).  URL https://link.aps.org/doi/10.1103/PhysRevLett.77.3787.
  * (50) Blanco-Redondo, A. _et al._ Observation of soliton compression in silicon photonic crystals.  _Nature Communications_ 5, 3160 (2014).  URL https://doi.org/10.1038/ncomms4160.
  * (51) Baboux, F. _et al._ Bosonic condensation and disorder-induced localization in a flat band.  _Phys. Rev. Lett._ 116, 066402 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevLett.116.066402.
  * (52) Amelio, I. & Carusotto, I.  Theory of the coherence of topological lasers.  _Phys. Rev. X_ 10, 041060 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevX.10.041060.
  * (53) Comaron, P., Shahnazaryan, V., Brzezicki, W., Hyart, T. & Matuszewski, M.  Non-hermitian topological end-mode lasing in polariton systems.  _Phys. Rev. Research_ 2, 022051 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.022051.
  * (54) Zamora, A., Sieberer, L. M., Dunnett, K., Diehl, S. & Szymańska, M. H.  Tuning across universalities with a driven open condensate.  _Phys. Rev. X_ 7, 041006 (2017).  URL https://link.aps.org/doi/10.1103/PhysRevX.7.041006.
  * (55) Fruchart, M., Hanai, R., Littlewood, P. B. & Vitelli, V.  Non-reciprocal phase transitions.  _Nature_ 592, 363–369 (2021).  URL https://doi.org/10.1038/s41586-021-03375-9.
  * (56) Mandal, S., Banerjee, R., Ostrovskaya, E. A. & Liew, T. C. H.  Nonreciprocal transport of exciton polaritons in a non-hermitian chain.  _Phys. Rev. Lett._ 125, 123902 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.123902.
  * (57) Kolkowski, R. & Koenderink, A. F.  Lattice resonances in optical metasurfaces with gain and loss.  _Proceedings of the IEEE_ 108, 795–818 (2020).  URL https://ieeexplore.ieee.org/abstract/document/8852732.

## Acknowledgements

We acknowledge the support of the UK’s Engineering and Physical Sciences
Research Council (grant EP/M025330/1 on Hybrid Polaritonics), the RFBR
projects No. 20-52-12026 (jointly with DFG) and No. 20-02-00919. H.S. and
P.G.L acknowledge support by the European Union’s Horizon 2020 program,
through a FET Open research and innovation action under the grant agreement
No. 899141 (PoLLoC).

## Author contributions

P.G.L. led the research project. P.G.L., S.A., A.A., J.D.T. designed the
experiment. S.A., A.A. and J.D.T. carried out the experiment and analysed the
data. H.S. developed the theoretical modelling and performed numerical
simulations. All authors contributed to the writing of the manuscript.

## Competing interests

The authors declare no competing interests.

Supplementary Information

## Supplementary Note 1: Additional experimental measurements on optical
lattices with square, Lieb, and inverse Lieb geometry

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\ \upmu$m (see
Supplementary Fig. S1) and Lieb lattice $D=22.3\ \upmu$m (see Supplementary
Fig. S2). We also present complementary data on the inverse Lieb lattice in
Fig. 4 in the main manuscript (see Supplementary Fig. S3). The full-width-
half-maximum (FWHM) of the Gaussian laser spots exciting the condensates is
$\approx 2$ $\upmu$m. Just like in the main text, we define the real space
horizontal and vertical coordinates as $\mathbf{r}=(x,y)$, and in momentum
space $\mathbf{k}=(k_{x},k_{y})$.

Supplementary Figs. S1a and S1b show the PL in real space and energy-momentum
space (along $k_{x}=0$) approximately at threshold ($P\approx P_{\text{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
Supplementary Figs. S1c and S1d we show the same measurements above
condensation threshold ($P\approx 1.3P_{\text{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
Supplementary 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.

Figure S1: Optical lattice of 25 polariton condensates arranged in a square
geometry with a lattice constant set to $D=22.5$ $\upmu$m. (a,c) Real space
and (b,d) energy-resolved momentum space polariton PL at ($P\approx
P_{\text{thr}}$) and above ($P\approx 1.3P_{\text{thr}}$) threshold
respectively.

In Supplementary Fig. S2 we provide similar measurements to those in Fig. 2a-e
in the main manuscript but using a larger lattice constant of $D=22.3\
\upmu$m. The results presented in Supplementary Fig. S2 show that by changing
the lattice constant in the 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 (main manuscript) 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.

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$ $\upmu$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\approx 1.1P_{\text{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\approx 1.1P_{\text{thr}}$. Semi-
transparent thin lines denoted with yellow letters “d” and “e” correspond to
$k_{x}=0$ and $k_{x}=2\pi/D$, for which the energy-resolved momentum space PL
shown in (d),(e) respectively.

In Supplementary Fig. S3 we provide experimental results for the inverse
(conventional tight binding) Lieb lattice hosting $P$-flatbands corresponding
to Figs. 4e-g in the main text. Supplementary Fig. S3a shows the energy-
resolved reciprocal space PL taken at $k_{x}=0$ where yellow arrows indicate
the location of $P$\- and $S$-bands. Supplementary Fig. S3b is measured at
$k_{x}=2\pi/D$ where the two horizontally dashed lines, denoted with yellow
letters “(c)” and “(d)” correspond to the momentum space isoenergy planes PL
shown in Supplementary Figs. S3c and S3d. We note that in each panel the PL is
independently normalised. The data shown in Supplementary 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
$\Gamma\leftrightarrow X$ edges in the outer Brillouin zones. The lowest
$P$-band is completely flat along the $\Gamma\leftrightarrow X$ edges but
possesses strong dispersion along the $X\leftrightarrow
M\leftrightarrow\Gamma$ edges as shown in Fig. 4i in the main manuscript. The
same type of PL is also observed in our calculations (see Supplementary Fig.
S7e).

Figure S3: Experimentally measured band structure in an inverse Lieb lattice
with $D=16$ $\upmu$m. (a,b) Energy-resolved momentum space polariton PL at
$k_{x}=0$ and $k_{x}=2\pi/D$ corresponding to the system shown in Figs. 4e-g
of the main text. Yellow arrows in (a) indicate $P$-band and $S$-band states
respectively. White dashed lines in (b) with yellow letters “(c)” and “(d)”
specify the energies of the momentum space PL presented in (c) and (d),
respectively. All panels are normalised independently.

## Supplementary Note 2: 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.

$i\hbar\frac{d\Psi}{dt}=\left[\hat{E}+V(\mathbf{r})\right]\Psi$ (S1)

Here, $\hat{E}$ 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 SM_kavokin2017microcavities ,

$\hat{E}(\mathbf{k})=\frac{1}{2}\left[\hat{E}_{\phi}(\mathbf{k})-\sqrt{\hat{E}_{\phi}^{2}(\mathbf{k})+(\hbar\Omega)^{2}}\right],$
(S2)

where $\hat{E}_{\phi}(\mathbf{k})$ is the cavity photon dispersion and
$\hbar\Omega=12$ meV is the Rabi energy of our sample
SM_cilibrizzi_polariton_2014 . The longitudinally confined photons obey the
dispersion relation,

$\hat{E}_{\phi}(\mathbf{k})=\sqrt{\left(\frac{hc}{\lambda}\right)^{2}+\left(\frac{\hbar
ck}{n}\right)^{2}}+\Delta-\frac{i\hbar\gamma}{2},$ (S3)

where $\lambda=860$ nm is the cavity resonance, $n=3.556$ is the cavity
refractive index, $\Delta=-5$ meV is the cavity photon detuning from the
exciton resonance, and $\gamma=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 Supplementary Fig.
S4a. In Supplementary Fig. S4b we show a schematic of the Lieb lattice and the
three square sublattices denoted by their label $A,B,C$.

Figure S4: Bare lower polariton dispersion and Lieb lattice illustration. (a)
Heatmap shows the normalized measured lower polariton dispersion with Eq. (S2)
plotted on top as a white dashed line. (b) Schematic of the Lieb lattice
geometry consisting of three square sublattices labeled $A,B,C$.

The potential $V(\mathbf{r})$ is directly proportional to the incident non-
resonant laser intensity and can be written as a superposition of spatially
separated Gaussians,

$V(\mathbf{r})=V_{0}\sum_{n=1}^{N}\exp{\left[-\frac{(x-x_{n})^{2}+(y-y_{n})^{2}}{2w^{2}}\right]}.$
(S4)

Here, $N$ is the number of pumps in question and $x_{n},y_{n}$ are their
coordinates, $V_{0}$ denotes the complex potential strength.
$\text{Re}{(V_{0})}>0$ and $\text{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$
$\upmu$m is the RMS width of the pump spots corresponding to 3 $\upmu$m FWHM.
The width is slightly larger than that of the lasers (2 $\upmu$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 $\Psi(\mathbf{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 $\tilde{\psi}(\omega,\mathbf{k})$ survive
the longest in the system (i.e., decay the slowest) in presence of optical
gain at the pump spots [i.e., $\text{Im}{(V_{0})}>0$]. Results of this method
are shown in Supplementary Figs. S5-S7 where we extract numerically the
dispersion from the dynamics of many random initial realizations of
$\tilde{\psi}{(\omega,\mathbf{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.

Figure S5: Calculated band structure in a square arranged scatterer lattice of
Gaussian spots. (a) Non-resonant laser profile $V(\mathbf{r})$ (normalised)
shaped into a square scatterer lattice of Gaussians with a lattice constant
$D=22.5$ $\upmu$m and $V_{0}=2+0.25i$ meV. (b,c) Energy-momentum space PL at
(b) $k_{y}=0$ and (c) $k_{y}=2\pi/D$. (d-f) Polariton momentum space PL around
threshold at $E=1.72,1.3,0.87$ meV respectively corresponding to the white
dotted lines in (b,c). Figure S6: Calculated band structure in a Lieb arranged
scatterer lattice of Gaussian spots. (a) Non-resonant laser profile
$V(\mathbf{r})$ (normalised) shaped into a Lieb scatterer lattice of Gaussians
with a lattice constant $D=20.3$ $\upmu$m and $V_{0}=2.2+0.25i$ meV
corresponding to Fig. 2 in the main text. (b,c) Energy-momentum space PL
around threshold for (b) $k_{y}=0$ and (c) $k_{y}=2\pi/D$. (d-f) Polariton
momentum space PL around threshold at $E=1.54,1.38,0.79$ meV respectively
corresponding to the dotted lines in (b,c). Figure S7: Calculated band
structure in a Lieb lattice formed from inverse arrangement of Gaussian pump
spots. (a) Non-resonant laser profile shaped to form a conventional Lieb
lattice of potential minima between the lasers. Lattice constant is $D=15.6$
$\upmu$m and $V_{0}=2$ meV corresponding to Fig. 4 in the main text. (b,c)
Energy-momentum space PL around threshold for (b) $k_{y}=0$ and (c)
$k_{y}=2\pi/D$. (d-f) Polariton momentum space PL around threshold at
$E=0.47,0.38,0.2$ meV respectively corresponding to the dotted lines in (b,c).

## Supplementary Note 3: 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,

$\hat{E}(\mathbf{k})=\frac{\hbar^{2}k^{2}}{2m}$ (S5)

where the effective mass of the lower polaritons is $m=0.32$ meV ps2
$\upmu$m-2 found by fitting the dispersion in Supplementary Fig. S4a. The
potential is assumed to be infinite and periodic such that
$V(\mathbf{r})=V(\mathbf{r}+\mathbf{D})$ where
$\mathbf{D}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}$ is the translational
symmetry vector defined in the bases of primitive lattice vectors
$\mathbf{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 $\mathbf{q}=(q_{x},q_{y})$ and Bloch
states in the $n$th band $u_{n,\mathbf{q}}(\mathbf{r})$,

$\Psi_{n}(\mathbf{r})=e^{i\mathbf{q}\cdot\mathbf{r}}u_{n,\mathbf{q}}(\mathbf{r}).$
(S6)

Here, $u_{n,\mathbf{q}}(\mathbf{r})=u_{n,\mathbf{q}}(\mathbf{r}+\mathbf{D})$.
Substituting Eq. (S6) into Eq. (S1) the Bloch eigenvalue problem then reads,

$\left[\frac{\hbar^{2}}{2m}\left[\left(q_{x}-i\frac{\partial}{\partial
x}\right)^{2}+\left(q_{y}-i\frac{\partial}{\partial
y}\right)^{2}\right]+V(\mathbf{r})\right]u_{n,\mathbf{q}}(\mathbf{r})=E_{n,\mathbf{q}}u_{n,\mathbf{q}}(\mathbf{r}).$
(S7)

We can diagonalise Eq. (S7) by expressing the Gaussian potential as a Fourier
series in the basis of the reciprocal lattice vectors defined through
$\mathbf{G}\cdot\mathbf{D}=2\pi n$ for $n\in\mathbb{Z}$,

$V(\mathbf{r})=\sum_{G}V_{G}e^{i\mathbf{G}\cdot\mathbf{r}}.$ (S8)

Supplementary Fig. S8a shows a portion of the calculated band structure for a
Lieb arrangement of the Gaussian spots (see Supplementary Fig. S8f) with a
lattice constant of $D=22.3$ $\upmu$m with some features of interest marked
with magenta colored boxes. Supplementary Fig. S8b shows the real space
intensity of the Bloch state $|u_{n_{c},\mathbf{q}=0}|^{2}$, marked by red
circle corresponding to the experimentally observed pattern shown in
Supplementary Fig. S2. A portion of the Gaussian pump locations are marked by
red squares in both Supplementary Figs. S8b and S8f for clarity. In
Supplementary Fig. S8a we observe regions in the band structure with almost
linear- and flat-like bands (magenta squares zoomed in Supplementary 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 Supplementary Fig.
S8b shows complete destructive interference at sublattice $B$ whereas in the
experimentally measured spatial PL in Supplementary 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 Supplementary 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).

Figure S8: Calculated polariton bands using Bloch’s theorem in the Lieb
scatterer lattice. (a) Calculated bands using Eq. (S7) for a Lieb lattice of
Gaussians, similar to Fig. S6 and Fig. S2. (b) Calculated Bloch state
intensity $|u_{n_{c},\mathbf{q}=0}(\mathbf{r})|^{2}$ belonging to the band
marked by the red circle in (a) and the condensate real space PL from Fig.
2(b) in the main text. The arrangement of the Gaussian pumps is indicated by
the red squares. (c-e) Zoomed in regions of the dispersion showing exotic
points where the dispersion becomes practically linear (massless polaritons)
or flat (heavy polaritons). (f) Potential profile $V(\mathbf{r})$ (normalised)
shaped into a Lieb scatterer lattice of Gaussians with a lattice constant
$D=22.3$ $\upmu$m and $V_{0}=1.25$ meV.

## Supplementary Note 4: 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 SM_Wouters_PRL2007  written for the
condensate wavefunction, $\Psi({\bf r},t),$ and a rate equation for the
exciton reservoir supplying particles into the condensate, $n({\bf r},t)$:

$\displaystyle i\frac{\partial\Psi}{\partial t}$
$\displaystyle=\Big{[}\frac{-\hbar\nabla^{2}}{2m}+\alpha|\Psi|^{2}-\frac{i\gamma}{2}+\left(g+\frac{iR}{2}\right)n+GP(\mathbf{r})\Big{]}\Psi+i\mathfrak{R}(\Psi),$
(S9) $\displaystyle\frac{\partial n}{\partial t}$
$\displaystyle=-(\Gamma_{R}+R|\Psi|^{2})n+P(\mathbf{r}).$ (S10)

Here, $m$ is the polariton effective mass, $\alpha$ is the polariton-polariton
interaction strength, $R$ is the scattering rate of reservoir excitons into
the condensate, $\gamma$ is the rate of polariton losses through the cavity
mirrors, $\Gamma_{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(\mathbf{r})$ (proportional to the laser
profile and power) that do not scatter into the condensate, respectively. The
final term $\mathfrak{R}(\Psi)$ is an energy relaxation term
SM_Wouters_PRB2010 ; SM_PhysRevB.92.035305  that assists ground state
condensation and is taken proportional to the background density of bottleneck
excitons,

$\mathfrak{R}(\Psi)=(\lambda
n)^{\xi}\left[\frac{\hbar\nabla^{2}}{2m}+\mu(\mathbf{r},t)\right]\Psi.$ (S11)

Here, $\lambda$ denotes the energy relaxation efficiency, $\xi$ is a fitting
parameter, and $\mu(\mathbf{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(\mathbf{r},t)$ away from the laser pump spots we approximate the reservoir
driving profile with a Lorentzian,

$P(\mathbf{r})=P_{0}\sum_{n=1}^{N}\frac{L^{2}}{L^{2}+(x-x_{n})^{2}+(y-y_{n})^{2}}.$
(S12)

Here, $N$ is the number of pumps in question and $x_{n},y_{n}$ are their
coordinates, $2L=2$ $\mu$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. Supplementary
Note 2: Theory of polariton dynamics 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) SM_Tosi_NatPhys2012  and stochastic noise
SM_toepfer_coherence2020 . The lower polariton mass and lifetime are taken
corresponding to the properties of our cavity: $m=0.32$ meV ps2 $\mu$m-2, and
$\gamma=1/5.5$ ps-1. The remaining parameters are taken similar to those used
in previous works: $\hbar\alpha=3.3$ $\mu$eV $\mu$m2, $g=\alpha$,
$G=2g/\Gamma_{R}$, $R=0.8\alpha$, $\Gamma_{R}=2\gamma$, $\xi=1/4$, and
$\lambda^{\xi}=0.004$ $\mu$m2.

In each simulation we apply stochastic initial conditions to
$\Psi(\mathbf{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
$\langle|\Psi(\mathbf{r},t)|^{2}\rangle$ 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
$|\Psi(\mathbf{r})|^{2}=|u_{3,\mathbf{q}_{X}}(\mathbf{r})|^{2}+|u_{3,\mathbf{q}_{X^{\prime}}}(\mathbf{r})|^{2}$
with $\mathbf{q}_{X}=(0,\pi/D)$ and $\mathbf{q}_{X^{\prime}}=(\pi/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.

## Supplementary Note 5: 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 $\approx 5$ $\upmu$eV different
in energy from both the $\Gamma$ or the $X$ points (see Supplementary Fig.
S9b). Interestingly, the curvature of one flatband corresponds to lattice
particles of negative effective mass and the other to positive effective mass.
Intuitively, the nonequilibrium nature and energy relaxation of polaritons
triggers condensation into the minima of the lower flatband. The source of
this curvature is attributed to the finite potential depth and the Gaussian
shape of the excitation beams (see Supplementary 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 Supplementary Fig. S9c. The flatbands
in Supplementary 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.

Figure S9: Effects of the inverse optical Lieb lattice potential height on the
P-bands. (a) Four unit cells of the Pump induced potential landscape used in
Bloch’s analysis. (b) Calculated band structure for peak potential height of 2
meV. (c) Calculated band structure for peak potential height of 4 meV. The
polariton mass is set to $m=0.32$ meV ps2 $\mu$m-2.

## Supplementary Note 6: 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,

$\displaystyle\hat{H}=$
$\displaystyle\sum_{n,m}\bigg{[}(\hat{b}^{(x)}_{n,m})^{\dagger}[J_{2}(\hat{a}_{n,m}^{(x)}+\hat{a}_{n-1,m}^{(x)})-J_{1}(\hat{c}_{n,m}^{(x)}+\hat{c}_{n,m-1}^{(x)})]+\text{h.c.}$
$\displaystyle+(\hat{b}^{(y)}_{n,m})^{\dagger}[J_{1}(\hat{a}_{n,m}^{(y)}+\hat{a}_{n-1,m}^{(y)})-J_{2}(\hat{c}_{n,m}^{(y)}+\hat{c}_{n,m-1}^{(y)})]+\text{h.c.}$
$\displaystyle+\epsilon_{a,x}(\hat{a}_{n,m}^{(x)})^{\dagger}\hat{a}_{n,m}^{(x)}+\epsilon_{b,x}(\hat{b}_{n,m}^{(x)})^{\dagger}\hat{b}_{n,m}^{(x)}+\epsilon_{c,x}(\hat{c}_{n,m}^{(x)})^{\dagger}\hat{c}_{n,m}^{(x)}$
$\displaystyle+\epsilon_{a,y}(\hat{a}_{n,m}^{(y)})^{\dagger}\hat{a}_{n,m}^{(y)}+\epsilon_{b,y}(\hat{b}_{n,m}^{(y)})^{\dagger}\hat{b}_{n,m}^{(y)}+\epsilon_{c,y}(\hat{c}_{n,m}^{(y)})^{\dagger}\hat{c}_{n,m}^{(y)}\bigg{]}.$
(S13)

Here, $\hat{a}^{(x,y)}_{n,m},\hat{b}^{(x,y)}_{n,m},\hat{c}^{(x,y)}_{n,m}$
denote the annihilation operators at sublattices $A,B,C$ respectively [see
Supplementary Fig. S4], $(n,m)$ denotes the site index of the corresponding
sublattice, and $(x,y)$ P-orbital lobe alignment along the two orthogonal
coordinates $x,y$ of the lattice. The operators obey bosonic commutation
relation
$[\hat{a}_{n,m}^{(\alpha)},(\hat{a}_{n^{\prime},m^{\prime}}^{(\beta)})^{\dagger}]=\delta_{\alpha,\beta}\delta_{n,n^{\prime}}\delta_{m,m^{\prime}}$
and $[\hat{a}_{n,m}^{(\alpha)},\hat{a}_{n^{\prime},m^{\prime}}^{(\beta)}]=0$.
The coefficients $J_{1,2}\in\mathbb{C}$ denote the polariton non-Hermitian
hopping amplitudes between P-orbitals with orbital axes parallel
($\sigma$-bond) and perpendicular ($\pi$-bond) to the bond axis respectively.
The coefficients $\epsilon_{N,\alpha}$ denote on-site energy of the two
distinct P-orbitals at each sublattice.

Performing the standard Fourier transform of the single particle operators one
can transform the problem from the basis of lattice sites ($n,m$) into the
basis of crystal momentum $\mathbf{q}=(q_{x},q_{y})$.

$\frac{\hat{H}_{\mathbf{q}}}{2}=\begin{pmatrix}\hat{H}_{P_{x}}&0_{3,3}\\\
0_{3,3}&\hat{H}_{P_{y}}\end{pmatrix},$ (S14)

where

$\displaystyle\hat{H}_{P_{x}}$
$\displaystyle=\begin{pmatrix}\epsilon_{a,x}&J_{2}e^{i\frac{q_{y}}{2}}\cos{(\frac{q_{y}}{2})}&0\\\
J_{2}e^{-i\frac{q_{y}}{2}}\cos{(\frac{q_{y}}{2})}&\epsilon_{b,x}&J_{1}e^{i\frac{q_{x}}{2}}\cos{(\frac{q_{x}}{2})}\\\
0&J_{1}e^{-i\frac{q_{x}}{2}}\cos{(\frac{q_{x}}{2})}&\epsilon_{c,x}\end{pmatrix},$
$\displaystyle\hat{H}_{P_{y}}$
$\displaystyle=\begin{pmatrix}\epsilon_{a,y}&J_{1}e^{i\frac{q_{y}}{2}}\cos{(\frac{q_{y}}{2})}&0\\\
J_{1}e^{-i\frac{q_{y}}{2}}\cos{(\frac{q_{y}}{2})}&\epsilon_{b,y}&J_{2}e^{i\frac{q_{x}}{2}}\cos{(\frac{q_{x}}{2})}\\\
0&J_{2}e^{-i\frac{q_{x}}{2}}\cos{(\frac{q_{x}}{2})}&\epsilon_{c,y}.\end{pmatrix}.$

Setting all diagonal terms to zero corresponds to all sites being in resonance
and with equal on-site gain. In Supplementary 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 $\Gamma$ point in the upper bands.

Figure S10: Calculated P-orbital bands using tight binding theorem. Blue to
red colorscale shows energies of small and large imaginary part (gain)
respectively. (a) All diagonal elements $\epsilon_{N,\alpha}=0$ resulting in
highest gain at the top of the conduction $P$-bands. (b) Non-zero diagonal
elements introducing splitting between $P_{x}$ and $P_{y}$ orbitals at
sublattices $A$ and $C$ results in gain moving towards the flatband.

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 Supplementary 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
$\epsilon_{a,x}=\epsilon_{c,y}=-\epsilon_{0}+i\kappa$,
$\epsilon_{a,y}=\epsilon_{c,x}=\epsilon_{0}-i\kappa$, and
$\epsilon_{b,x}=\epsilon_{b,y}=0$ where $\epsilon_{0},\kappa>0$. Such
splitting arises from the slightly elliptically shaped $A$ and $C$ sublattice
sites [see Supplementary Fig. S7a]. Physically, $\kappa$ captures the
increased gain of the lower energy state due to greater overlap with the gain
region from the Gaussian pumps. In contrast to Supplementary Fig. S10a the
flatband between the $X$ and the $M$ points in Supplementary 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.

## References

  * (1) Kavokin, A. V., Baumberg, J. J., Malpuech, G. & Laussy, F. P.  _Microcavities_ , vol. 21 (Oxford university press, 2017).
  * (2) Cilibrizzi, P. _et al._ Polariton condensation in a strain-compensated planar microcavity with ingaas quantum wells.  _Applied Physics Letters_ 105, 191118 (2014).  URL https://doi.org/10.1063/1.4901814.
  * (3) Wouters, M. & Carusotto, I.  Excitations in a nonequilibrium bose-einstein condensate of exciton polaritons.  _Phys. Rev. Lett._ 99, 140402 (2007).  URL https://link.aps.org/doi/10.1103/PhysRevLett.99.140402.
  * (4) Wouters, M., Liew, T. C. H. & Savona, V.  Energy relaxation in one-dimensional polariton condensates.  _Phys. Rev. B_ 82, 245315 (2010).  URL https://link.aps.org/doi/10.1103/PhysRevB.82.245315.
  * (5) Askitopoulos, A. _et al._ Robust platform for engineering pure-quantum-state transitions in polariton condensates.  _Phys. Rev. B_ 92, 035305 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevB.92.035305.
  * (6) Tosi, G. _et al._ Sculpting oscillators with light within a nonlinear quantum fluid.  _Nature Physics_ 8, 190–194 (2012).  URL https://doi.org/10.1038/nphys2182.
  * (7) Töpfer, J. D. _et al._ Engineering spatial coherence in lattices of polariton condensates.  _Optica_ 8, 106–113 (2021).  URL http://www.osapublishing.org/optica/abstract.cfm?URI=optica-8-1-106.


# XMM-Newton detection of soft time lags in the TDE candidate AT 2018fyk

Wenda Zhang
National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun
Road, Beijing 100101, China E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this letter we report a tentative detection of soft time lags (i.e.
variability of softer photons lags behind the variability of harder photons)
in one XMM-Newton observation of the tidal disruption event (TDE) candidate AT
2018fyk while the source was in the hard spectral state. The lags are detected
at $6.51\times 10^{-5}~{}\rm Hz$. The amplitude of the lags with respect to
0.5$-$1 keV monotonically decreases with the photon energy, from $\sim
1200~{}\rm s$ at 0.3$-$0.5 keV to $\sim-4200~{}\rm s$ at 3$-$5 keV (in our
convention a positive lag means lagging behind the reference band). We find
that the amplitude is proportional to the logarithm of the energy separation
between the examined band and the reference band. The energy-dependent
covariance spectrum indicates that the correlated variability is more likely
to be associated with the non-thermal radiation. The soft lags are difficult
to reconcile with the reverberation scenario that are used to explain the soft
lags in active galactic nuclei. On the other hand, the observed soft lags are
consistent with the picture that the soft X-rays are down-scattered hard
X-rays by the outflow as predicted by “unification” models of TDEs.

###### keywords:

transients: tidal disruption events – accretion, accretion discs – black hole
physics

††pagerange: XMM-Newton detection of soft time lags in the TDE candidate AT
2018fyk–Appendix

##  1 Introduction

Tidal disruption events (TDEs) are transient events caused by super-massive
black holes (SMBHs) tidally disrupting nearby stars and accreting stellar
debris (see Gezari, 2021, for a review). TDEs offer us unique opportunities to
study dormant SMBHs. So far, a few tens of TDEs have been discovered.
Depending on their X-ray properties, TDEs can be classified into two
categories: X-ray TDEs that are detected in X-rays, and optical-UV TDEs that
lack X-ray emission. Understanding the X-ray/optical-UV dichotomy is an
important task in the study of TDEs.

To explain the weak X-ray emission of optical-UV TDEs, several models have
been proposed. The first set of models (often referred to as “unification”
models of TDEs) claimed that the dichotomy is a viewing angle effect: the
X-ray and optical-UV TDEs are intrinsically the same, but for observers at
large inclinations the X-ray radiation is obscured by the wind launched by the
super-Eddington accretion disc (e.g. Dai et al., 2018). In the second set of
models, the dichotomy is intrinsic; the strength of the X-ray emission depends
on the physical properties of the system. For instance, it has been proposed
that the disc formation is faster for TDEs with larger black hole mass (e.g.
Guillochon & Ramirez-Ruiz, 2015), resulting in the absence of X-ray emission
in the early outburst stage of low-mass TDEs.

AT 2018fyk is a candidate X-ray TDE discovered by the ASAS-SN on September
8th, 2018 (Stanek, 2018). It is located at the nucleus of LCRS
B224721.6$-$450748, a galaxy at a redshift of
0.06111https://ned.ipac.caltech.edu/byname?objname=LCRS+B224721.6-450748.
Wevers (2020) measured the mass of the central SMBH of the galaxy using the
$M-\sigma$ relation, and found ${\rm log_{10}}\ M_{\rm BH}=7.7\pm 0.4$ where
$M_{\rm BH}$ is the black hole mass in solar mass. Wevers et al. (2019)
performed a multi-wavelength study of AT 2018fyk and found evidences for rapid
disc formation. Wevers et al. (2021) studied the long-term evolution of the
X-ray emission of AT 2018fyk, and discovered an X-ray state transition of AT
2018fyk from a soft, thermal-dominated spectral state to a hard, non-thermal
dominated state. The spectral state transition is accompanied by a change in
the timing properties, with high-frequency variability ($\gtrsim 10^{-5}~{}\rm
Hz$) emerging after the state transition.

We perform Fourier time-lag analysis of AT 2018fyk. Fourier time-lag analysis
has been proved to be a great tool to study the accretion flow around black
holes with mass across a large dynamic range. Time-lag analysis of both active
galactic nuclei (AGNs) and black hole X-ray binaries (BHXRBs) has revealed
continuum hard lags (e.g. Papadakis et al., 2001; McHardy et al., 2004;
Arévalo et al., 2006; Arévalo et al., 2008; Sriram et al., 2009; Epitropakis &
Papadakis, 2017; Papadakis et al., 2019; Miyamoto et al., 1988; Miyamoto &
Kitamoto, 1989; Nowak & Vaughan, 1996; Nowak et al., 1999; Cui et al., 2000;
Crary et al., 1998; Reig et al., 2003; Altamirano & Méndez, 2015). On top of
the continuum hard lags, there have been detections of soft lags where the
variability of softer bands are lagging behind the harder bands. The soft lags
were first detected by Fabian et al. (2009) in 1H 0707$-$495, and then in
other AGNs (e.g. Zoghbi & Fabian, 2011; Cackett et al., 2013; Fabian et al.,
2013; De Marco et al., 2013a; Alston et al., 2014), as well as in BHXRBs (e.g.
Uttley et al., 2011; De Marco et al., 2015; De Marco et al., 2017). The soft
lags are believed to be the lags between the hard coronal emission and the
reflection emission due to the hard coronal emission irradiating the disc, and
are consequently referred to as “reverberation lags”. De Marco et al. (2013a)
performed a systematic study of the reverberation lags, and found correlations
between the amplitude of the reverberation lags and the black hole mass, and
between the variability frequency and the black hole mass.

Soft lags have been detected in ultra-luminous X-ray sources (ULXs) as well
(Heil & Vaughan, 2010; De Marco et al., 2013b; Hernández-García et al., 2015;
Pinto et al., 2017; Kara et al., 2020; Pintore et al., 2021; Mondal et al.,
2021). ULXs are extragalactic, off-nuclear, ultra-luminous ($\gtrsim
10^{39}~{}\rm erg~{}s^{-1}$) X-ray sources. It is generally believed that the
majority of ULXs are stellar-mass black holes or neutron stars accreting at
super-Eddington rates, instead of intermediate-mass black holes (IMBHs)
accreting at sub-Eddington rates. Most of the soft lags in ULXs have
amplitudes of a few hundreds to thousands of seconds and are usually
interpreted as the time lags due to hard X-ray photons getting down-scattered
by the outflow.

In this letter we report a detection of soft lags in AT 2018fyk. The letter is
organised as follows. The procedure of data reduction is illustrated in Sec.
2. The results are presented in Sec. 3. We discuss the results in Sec. 4.

##  2 Data Reduction

XMM-Newton observed AT 2018fyk three times. The source was detected in the
first (obs id: 0831790201) and second (obs id: 0853980201) observations, and
undetected in the last observation. During the first observation AT 2018fyk
was in the soft state, and no X-ray variability was detected above $\sim
10^{-5}~{}\rm Hz$, while during the second observation strong high-frequency
($\gtrsim 10^{-5}~{}\rm Hz$) X-ray variability was discovered (Wevers et al.,
2021). We therefore perform Fourier time-lag analysis on the pn data of
0853980201.

The observation 0853980201 were taken on Oct. 27th, 2019. During the
observation pn was operating in the full-frame mode. The data are reduced with
the the XMM-Newton science analysis system (SAS) version 18 following standard
procedures222https://www.cosmos.esa.int/web/xmm-newton/sas-threads. We produce
the pn clean event list with epproc. We filter the event list to remove
periods contaminated by high particle background flares, where the good time
interval is defined to be the period during which the 10$-$12 keV count rate
is no larger than $0.4~{}\rm counts/s$333https://www.cosmos.esa.int/web/xmm-
newton/sas-thread-epic-filterbackground. We are left with a $\sim 30~{}\rm ks$
exposure after removing flaring backgrounds. The source and background
lightcurves are extracted from the filtered event list with evselect. When
extracting the lightcurves we select single and double events ($PATTERN\leq
4$) with $FLAG==0$. The source region is a circle centering on the source
position with a radius of 35”. The background lightcurves are extracted from a
source-free region close to the source with the same radius. The source
lightcurves are then corrected for background and various other effects using
epiclccorr. We also extract source and background spectra with evselect, and
generate corresponding ancillary file and redistribution matrix using arfgen
and rmfgen, respectively.

We extract 30-s lightcurves of 5 energy bands: $0.3-0.5~{}\rm keV$,
$0.5-1~{}\rm keV$, $1-2~{}\rm keV$, $2-3~{}\rm keV$, and $3-5~{}\rm keV$.
Above $\sim 5$ keV the background is stronger than the source so we exclude
the $>5~{}\rm keV$ bands in our analysis. For each energy band we split the
lightcurve into 2 segments, each $\sim 15~{}\rm ks$ long. Taking $0.5-1~{}\rm
keV$ as the reference, we follow Epitropakis & Papadakis (2016) to compute the
frequency-dependent Fourier time lags. We estimate the significance and
uncertainty of the time-lag measurements by performing Monte$-$Carlo
simulations. Detailed procedures of the simulations can be found in the
Appendix.

##  3 Results

###  3.1 Time-lag analysis

The frequency-resolved time lags and the coherence are presented in the upper
and lower panels of Fig. 1, respectively. Note that in our convention a
positive lag means lagging behind the reference band. The coherence is close
to unity at the lowest frequency $6.51\times 10^{-5}~{}\rm Hz$, decreases with
the frequency, and drops to zero above $\sim 2\times 10^{-3}~{}\rm Hz$.

The time lags are close to zero except at $6.51\times 10^{-5}~{}\rm Hz$, the
lowest frequency we can reach. At $6.51\times 10^{-5}~{}\rm Hz$, the
confidence levels of the time-lags are $92.0\%$, $97.5\%$, $97.9\%$, and
$97.9\%$, for 0.3$-$0.5 keV, 1$-$2 keV, 2$-$3 keV, and 3$-$5 keV,
respectively. The maximum confidence level of $97.9\%$ corresponds to a
significance of $\sim 1.8\sigma$. The most surprising result is that the time
lags are soft, i.e. variability of soft photons are lagging behind variability
of hard photons. This is opposite to AGNs where the continuum time lags are
hard (e.g. Arévalo et al., 2006; Arévalo et al., 2008; Sriram et al., 2009;
Epitropakis & Papadakis, 2017). The amplitude of the time lags monotonically
decreases with the energy, from $1200~{}\rm s$ at the 0.3$-$0.5 keV band to
$-4200~{}\rm s$ at the 3$-$5 keV band.

In Fig. 2 we plot the amplitude of the time lags at $6.51\times 10^{-5}~{}\rm
Hz$ as a function of the centroid energy of the energy bands. The data points
seem to fall on a straight line on the plot where the time lags and energy are
plotted in linear and logarithmic scales, respectively. This implies that the
time lags depend logarithmically on the energy separation between two energy
bands. We fit the data with a logarithmic model, and find:

$\tau=(-5.63\pm 0.80)\times 10^{3}\ {\rm
log_{10}}\left(\frac{E_{2}}{E_{1}}\right)~{}\rm s,$ (1)

where $\tau$ is the amplitude of the time lag, and $E_{1}$, $E_{2}$ are the
centroid energies of the two energy bands. The best-fit model is plotted in a
dashed line in Fig. 2. By comparing the data and the model, it is obvious that
the model fits the data well, confirming the logarithmic dependence.

Figure 1: The Fourier time lags (upper panel) with respect to $0.5-1~{}\rm
keV$ and the coherence (lower panel). The results of different energy bands
are plotted in different colors, as indicated in the plot. In our convention a
positive lag means lagging behind the reference band. For $3-5~{}\rm keV$, the
90% and 95% confidence levels are plotted in dashed and dotted lines,
respectively.  Figure 2: The time lags at $6.51\times 10^{-5}~{}\rm Hz$ as a
function of the central energy of the energy bands. The error bars correspond
to 68% uncertainty. We fit the lag-energy data with a logarithmic function and
the best-fit model is shown by the black dashed line.

###  3.2 Energy and covariance spectrum

We fit the pn energy spectrum with Xspec 12.11.1. The spectrum is rebinned
with specgroup to ensure minimum counts of 30 per bin. We fit the model to the
data using the $\chi^{2}$ statistics. We fit the pn spectrum with the
following spectral model: $phabs*(blackbody+powerlaw)$, which includes two
components: a thermal blackbody component, and a nonthermal powerlaw
component, both of which are absorbed by interstellar media. The spectrum and
best-fit model are plotted in Fig. 3. The best-fit model (black solid) yields
$\chi^{2}_{\nu}=1.47$ (with 83 degrees of freedom), and consists of a
dominating powerlaw component (black dashed) with
$\Gamma=2.32^{+0.09}_{-0.07}$ and a soft thermal component (black dotted) with
a temperature of $100^{+7}_{-4}~{}\rm eV$.

To investigate the origin of the correlated variability, we also calculate the
covariance spectrum in the frequency range of $3.7\times 10^{-5}-1.7\times
10^{-4}~{}\rm Hz$. We extract lightcurves with time resolution of $3000~{}\rm
s$, and split the lightcurves into 2 segments, each of which contains 9 bins.
Then we compute the covariance spectrum and its uncertainty following
(Wilkinson & Uttley, 2009), taking the 0.5$-$1 keV band as the reference. The
covariance spectrum is plotted in blue in Fig. 3. The shape of the covariance
spectrum indicates that the correlated variability is more likely to be
associated with the non-thermal emission than the thermal emission.

Figure 3: The covariance spectrum in blue colors. We also plot the energy
spectra (black dots), the best-fit model (black solid), and the two components
of the model: blackbody (black dotted) and powerlaw (black dashed).

##  4 Discussion and Conclusions

We tentatively detect soft time lags in one XMM-Newton observation of the TDE
candidate AT 2018fyk. The observation (obs id: 0853980201) was taken while the
source was in the hard spectral state. The lags are detected at $6.51\times
10^{-5}~{}\rm Hz$, and the amplitude of the lags with respect to the 0.5$-$1
keV reference band are proportional to the logarithm of the energy separation
between the examined band and the reference band, ranging from $1200$ to
$-4200$ seconds. The time-lags of 2$-$3 keV and 3$-$5 keV are most
significant, both with a confidence level of $97.9\%$.

A model fit confirms the apparent logarithmic dependence. The coherence is
close to unity at $6.51\times 10^{-5}~{}\rm Hz$ and drops rapidly with the
variability frequency.

We extract the energy spectrum and fit it with a two-component model. The
best-fit model consists of a thermal component with a temperature of $\sim
100~{}\rm eV$ and a dominating non-thermal component with a photon index of
$\sim 2$. We also calculate the energy-dependent covariance spectrum. The
shape of the covariance spectrum indicates that the correlated variability is
associated with the non-thermal emission.

Below we discuss possible physical origins of the soft lags.

###  4.1 Reverberation lag

It is well established that the soft lags in AGNs are “reverberation lags”,
i.e. the time lags between the primary hard X-ray radiation emitted by hot
coronae, and the reflection radiation of the accretion disc due to primary
radiation irradiating the disc. De Marco et al. (2013a) performed a systematic
study of the reverberation lags in AGNs, and found correlations between the
amplitude of the reverberation lag and the mass of the super-massive black
hole, and between the variability frequency and the black hole mass. The
authors found that, between 0.3$-$1 keV and 1$-$5 keV,

$\displaystyle{\rm log_{10}}\ \nu_{\rm lag}$ $\displaystyle=$
$\displaystyle-3.50\pm 0.07-0.47\pm 0.09\ {\rm log_{10}}\ M_{7},$ (2)
$\displaystyle{\rm log_{10}}\ |\tau|$ $\displaystyle=$ $\displaystyle 1.98\pm
0.08+0.59\pm 0.11\ {\rm log_{10}}\ M_{7},$ (3)

where $\nu_{\rm lag}$ is the frequency of the reverberation lag in Hz, $\tau$
is the amplitude of the reverberation lag in second, and $M_{7}$ is the mass
of the super-massive black hole in $10^{7}~{}\rm M_{\odot}$. Substituting the
black hole mass of AT 2018fyk (${\rm log_{10}}~{}M_{\rm BH}=7.7\pm 0.4$;
Wevers, 2020), we obtain ${\rm log_{10}}\ \nu_{\rm lag}=-3.83\pm 0.21$, and
${\rm log_{10}}\ |\tau|=2.39\pm 0.26$. For comparison, the soft lags in AT
2018fyk are detected at $6.51\times 10^{-5}~{}\rm Hz{}$, and the amplitude of
the time lags between 0.3$-$1 and 1$-$5 keV is expected to be $-3.74\pm
0.53~{}\rm ks$ (using Eq. 1). Either the variability frequency or the
amplitude is inconsistent with the reverberation lag in AGNs. This indicates
that the soft lags in AT 2018fyk are not reverberation lags. The reverberation
lag scenario is also disfavored by the absence of soft excess in the
covariance spectrum (Uttley et al., 2011).

###  4.2 Scattering

Alternatively, the soft lags could be due to hard photons getting down-
scattered by the materials in the outflow. The unification models (e.g. Dai et
al., 2018) predicted that, at intermediate inclinations, the observers are
still expected to see X-ray emission that originates from the outflow down-
scattering hard X-ray photons. Naturally we expect to see soft lags as on
average the photons lose energy by down-scattering.

In the Compton down-scattering process, the average energy loss per scattering
is

$\frac{<x^{\prime}-x>}{x}=-x,\quad{\rm if}~{}kT_{e}/m_{e}c^{2}\ll x\ll 1,$ (4)

where $x^{\prime}\equiv E^{\prime}/m_{e}c^{2}$ and $x\equiv E/m_{e}c^{2}$ are
the dimensionless energies of the photons after and before the scattering,
respectively, $m_{e}$ is the electron rest mass, and $kT_{e}$ is the
temperature of the down-scattering medium. For a photon with an energy of
$x_{0}(\ll 1)$, after $N$ scattering, the expected energy of the photon is

$x_{N}\approx x_{0}-Nx_{0}^{2}.$ (5)

Therefore

$\frac{x_{N}}{x_{0}}\approx 1-Nx_{0}^{2}\approx e^{-Nx_{0}^{2}},$ (6)

leading to

$N\approx{\rm ln}\left(\frac{x_{0}}{x_{N}}\right)/x_{0}^{2}.$ (7)

Naturally we expect the amplitude of the time lags to be proportional to the
number of scatterings. Therefore in this down-scattering scenario the time
lags between two energy bands are expected to be proportional to the logarithm
of the energy separation between two bands. Thus this scenario is also
consistent with the observations that the amplitude of the lag depends
logarithmically on the energy separation between two bands (Fig. 2).

It is interesting to note that soft lags with amplitudes of a few hundreds to
thousands of seconds have been detected in four ULXs: NGC 55 ULX1 (Pinto et
al., 2017), NGC 1313 X1 (Kara et al., 2020), NGC 4559 X7 (Pintore et al.,
2021), and NGC 7456 ULX-1 (Mondal et al., 2021). For the lag amplitude to be
consistent with the reverberation lags in AGNs, the black hole masses have to
be at least a few times $10^{7}~{}\rm M_{\odot}$ (De Marco et al., 2013a),
much too massive for ULXs: actually most ULXs are believed to contain stellar-
mass black holes or neutron stars instead of IMBHs. These time lags have been
interpreted as due to hard photons getting down-scattered by outflows, the
same scenario with what we use to explain the soft lags in AT 2018fyk. As both
TDEs and ULXs are accreting at super-Eddington rates where launching of
outflow is expected, it is not surprising that we detect soft lags due to the
same physical process in these two kinds of objects.

## acknowledgements

The author is grateful to the anonymous referee for his/her very useful
suggestions. The author would like to thank Iossif Papadakis and Chichuan Jin
for helpful discussion. The author acknowledges the support by the Strategic
Pioneer Program on Space Science, Chinese Academy of Sciences through grant
XDA15052100. This research has made use of data obtained through the High
Energy Astrophysics Science Archive Research Center Online Service, provided
by the NASA/Goddard Space Flight Center. This research makes uses of
Matplotlib (Hunter, 2007), a Python 2D plotting library which produces
publication quality figures.

## Data Availability

The data underlying this article are publicly available from NASA’s High
Energy Astrophysics Science Archive Research Center (HEASARC) archive at
https://heasarc.gsfc.nasa.gov and ESA’s XMM-Newton Science Archive (XSA) at
https://www.cosmos.esa.int/web/xmm-newton/xsa.

## References

  * Alston et al. (2014) Alston W. N., Done C., Vaughan S., 2014, MNRAS, 439, 1548
  * Altamirano & Méndez (2015) Altamirano D., Méndez M., 2015, MNRAS, 449, 4027
  * Arévalo et al. (2006) Arévalo P., Papadakis I. E., Uttley P., McHardy I. M., Brinkmann W., 2006, MNRAS, 372, 401
  * Arévalo et al. (2008) Arévalo P., McHardy I. M., Summons D. P., 2008, MNRAS, 388, 211
  * Cackett et al. (2013) Cackett E. M., Fabian A. C., Zogbhi A., Kara E., Reynolds C., Uttley P., 2013, ApJ, 764, L9
  * Crary et al. (1998) Crary D. J., Finger M. H., Kouveliotou C., van der Hooft F., van der Klis M., Lewin W. H. G., van Paradijs J., 1998, ApJL, 493, L71
  * Cui et al. (2000) Cui W., Zhang S. N., Chen W., 2000, ApJL, 531, L45
  * Dai et al. (2018) Dai L., McKinney J. C., Roth N., Ramirez-Ruiz E., Miller M. C., 2018, ApJL, 859, L20
  * De Marco et al. (2013a) De Marco B., Ponti G., Cappi M., Dadina M., Uttley P., Cackett E. M., Fabian A. C., Miniutti G., 2013a, MNRAS, 431, 2441
  * De Marco et al. (2013b) De Marco B., Ponti G., Miniutti G., Belloni T., Cappi M., Dadina M., Muñoz-Darias T., 2013b, MNRAS, 436, 3782
  * De Marco et al. (2015) De Marco B., Ponti G., Muñoz-Darias T., Nandra K., 2015, ApJ, 814, 50
  * De Marco et al. (2017) De Marco B., et al., 2017, MNRAS, 471, 1475
  * Epitropakis & Papadakis (2016) Epitropakis A., Papadakis I. E., 2016, A&A, 591, A113
  * Epitropakis & Papadakis (2017) Epitropakis A., Papadakis I. E., 2017, MNRAS, 468, 3568
  * Fabian et al. (2009) Fabian A. C., et al., 2009, Nature, 459, 540
  * Fabian et al. (2013) Fabian A. C., et al., 2013, MNRAS, 429, 2917
  * Gezari (2021) Gezari S., 2021, arXiv e-prints, p. arXiv:2104.14580
  * Guillochon & Ramirez-Ruiz (2015) Guillochon J., Ramirez-Ruiz E., 2015, ApJ, 809, 166
  * Heil & Vaughan (2010) Heil L. M., Vaughan S., 2010, MNRAS, 405, L86
  * Hernández-García et al. (2015) Hernández-García L., Vaughan S., Roberts T. P., Middleton M., 2015, MNRAS, 453, 2877
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Kara et al. (2020) Kara E., et al., 2020, MNRAS, 491, 5172
  * McHardy et al. (2004) McHardy I. M., Papadakis I. E., Uttley P., Page M. J., Mason K. O., 2004, MNRAS, 348, 783
  * Miyamoto & Kitamoto (1989) Miyamoto S., Kitamoto S., 1989, Nature, 342, 773
  * Miyamoto et al. (1988) Miyamoto S., Kitamoto S., Mitsuda K., Dotani T., 1988, Nature, 336, 450
  * Mondal et al. (2021) Mondal S., Różańska A., De Marco B., Markowitz A., 2021, MNRAS, 505, L106
  * Nowak & Vaughan (1996) Nowak M. A., Vaughan B. A., 1996, MNRAS, 280, 227
  * Nowak et al. (1999) Nowak M. A., Vaughan B. A., Wilms J., Dove J. B., Begelman M. C., 1999, ApJ, 510, 874
  * Papadakis et al. (2001) Papadakis I. E., Nandra K., Kazanas D., 2001, ApJL, 554, L133
  * Papadakis et al. (2019) Papadakis I. E., Rigas A., Markowitz A., McHardy I. M., 2019, MNRAS, 485, 1454
  * Pinto et al. (2017) Pinto C., et al., 2017, MNRAS, 468, 2865
  * Pintore et al. (2021) Pintore F., et al., 2021, MNRAS, 504, 551
  * Reig et al. (2003) Reig P., Kylafis N. D., Giannios D., 2003, A&A, 403, L15
  * Sriram et al. (2009) Sriram K., Agrawal V. K., Rao A. R., 2009, ApJ, 700, 1042
  * Stanek (2018) Stanek K. Z., 2018, Transient Name Server Discovery Report, 2018–1325, 1
  * Timmer & Koenig (1995) Timmer J., Koenig M., 1995, A&A, 300, 707
  * Uttley et al. (2011) Uttley P., Wilkinson T., Cassatella P., Wilms J., Pottschmidt K., Hanke M., Böck M., 2011, MNRAS, 414, L60
  * Wevers (2020) Wevers T., 2020, MNRAS, 497, L1
  * Wevers et al. (2019) Wevers T., et al., 2019, MNRAS, 488, 4816
  * Wevers et al. (2021) Wevers T., et al., 2021, ApJ, 912, 151
  * Wilkinson & Uttley (2009) Wilkinson T., Uttley P., 2009, MNRAS, 397, 666
  * Zoghbi & Fabian (2011) Zoghbi A., Fabian A. C., 2011, MNRAS, 418, 2642

## Appendix

To estimate the significance and uncertainty of the measured time-lags, for
each energy band (out of 0.3$-$0.5 keV, 0.5$-$1 keV, 1$-$2 keV, 2$-$3 keV, and
3$-$5 keV), we fit the Leahy-normalized power spectral density with the
following model:

$P(\nu)=N\nu^{-\alpha}+2,$ (8)

i.e. a powerlaw function plus a constant power of 2. The best-fit parameters
and the covariance matrix are obtained by minimizing $\chi^{2}$. Then, we
perform 10,000 simulations. In each simulation, for every energy band we
sample a powerlaw power spectral density using the best-fit model parameters
and the covariance matrix, simulate the lightcurve following Timmer & Koenig
(1995), and add Poisson noise. The phase of the Fourier transform is kept to
the the same with that of the reference band to ensure zero time-lag between
different energy bands. For the simulations we take a time resolution of 30 s,
the same with the real lightcurves. To avoid red noise leakage, the duration
of the simulated lightcurves is $\sim 180~{}\rm ks$, 5 times longer than the
real lightcurves. We randomly select a $\sim 30~{}\rm ks$ interval within the
whole $\sim 180~{}\rm ks$ interval, and apply identical procedures as in Sec.
2 to the $30~{}\rm ks$ time series to compute the time-lags for the four
energy bands. Finally we derive the confidence levels and uncertainties from
the distribution of the simulated time-lags.