# Engineering spatial coherence in lattices of polariton condensates

J. D. Töpfer<EMAIL_ADDRESS>Skolkovo Institute of Science and
Technology, Moscow, Bolshoy Boulevard 30, bld. 1, 121205, Russia  School of
Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ,
United Kingdom  I. Chatzopoulos  School of Physics and Astronomy, University
of Southampton, Southampton, SO17 1BJ, United Kingdom  H. Sigurdsson  Skolkovo
Institute of Science and Technology, Moscow, Bolshoy Boulevard 30, bld. 1,
121205, Russia  School of Physics and Astronomy, University of Southampton,
Southampton, SO17 1BJ, United Kingdom  T. Cookson  School of Physics and
Astronomy, University of Southampton, Southampton, SO17 1BJ, United Kingdom
Y. G. Rubo  Instituto de Energías Renovables, Universidad Nacional Autónoma de
México, Temixco, Morelos, 62580, Mexico  P. G. Lagoudakis
<EMAIL_ADDRESS>Skolkovo Institute of Science and Technology, Moscow,
Bolshoy Boulevard 30, bld. 1, 121205, Russia  School of Physics and Astronomy,
University of Southampton, Southampton, SO17 1BJ, United Kingdom

(August 27, 2024)

Artificial lattices of coherently coupled macroscopic states are at the heart
of applications ranging from solving hard combinatorial optimisation problems
to simulating complex many-body physical systems. The size and complexity of
the problems scales with the extent of coherence across the lattice. Although
the fundamental limit of spatial coherence depends on the nature of the
couplings and lattice parameters, it is usually engineering constrains that
define the size of the system. Here, we engineer polariton condensate lattices
with active control on the spatial arrangement and condensate density that
result in near-diffraction limited emission, and spatial coherence that
exceeds by nearly two orders of magnitude the size of each individual
condensate. We utilise these advancements to unravel the dependence of spatial
correlations between polariton condensates on the lattice geometry.

##  I Introduction

Synchronisation and the emergence of coherence between coupled elements are
universal concepts arising in nature and technology Pikovsky  _et al._ (2003).
They dictate collective human behaviour Strogatz  _et al._ (2005), functioning
of neurological systems Gray  _et al._ (1989), as well as phase transitions of
quantum systems to macroscopically collective entities at low temperatures
Davis  _et al._ (1995). For small-size systems such as two mechanical
pendulums coupled through a common support Bennett  _et al._ (2002) or two
coupled laser cavities Soriano  _et al._ (2013) frequency locking of the two
elements depends on the inter-element coupling strength in competition with
any inherent dephasing mechanisms. In larger systems consisting of many
interacting elements, such as social structures or coupled laser networks the
underlying coupling topology (or network architecture) critically influences
the dynamics and coherence formed in these systems Watts and Strogatz (1998);
Strogatz (2001). Engineering spatial coherence in such large networks is a key
element for increasing system performances in power grids Motter  _et al._
(2013), novel computational devices for classification tasks Romera  _et al._
(2018) or laser arrays for creation and control of high power beams Kao  _et
al._ (2016).

Lattices of coupled condensates are investigated for the simulation and
computation of complex tasks in atomic Morsch and Oberthaler (2006); Struck
_et al._ (2011); Bloch  _et al._ (2012), photonic Dung  _et al._ (2017) and
polaritonic Amo and Bloch (2016); Berloff  _et al._ (2017); Ballarini  _et
al._ (2020) platforms. Functionality and computational performance of these
systems is ultimately limited by the system’s spatial coherence length, i.e.
how many condensates can coherently be coupled. From a more fundamental point
of view, increasing the coherence allows one to construct a larger many-body
system and approach the ideal limit of a homogeneous “infinite” order
parameter, which is essential for the study of phase transitions in
interacting bosonic systems such as the Berezinskii-Kosterlitz-Thouless
transition in two-dimensions (2D) Hadzibabic  _et al._ (2006). In the case of
exciton-polariton condensates in semiconductor microcavites Kasprzak  _et al._
(2006), measurements were initially limited to small condensate sizes Deng
_et al._ (2007); Krizhanovskii  _et al._ (2009); Manni  _et al._ (2011);
Roumpos  _et al._ (2012) but advancements in sample fabrication and
experimental techniques now allow studies on their coherence properties as
extended objects far away from the excitation area Wertz  _et al._ (2010);
Caputo  _et al._ (2018), in optical traps Ohadi  _et al._ (2018); Askitopoulos
_et al._ (2019), and photonic microstructures Baboux  _et al._ (2018). The
emergence of correlations beyond the spatial extension of the laser excitation
area has underpinned the creation of networks and lattices of polariton
condensates.

Here, we engineer control over the particle density and position of each
polariton condensate across a polariton lattice under non-resonant optical
excitation. We overcome the challenge of disorder induced localisation and
dephasing by employing a feedback scheme for each individual laser spot,
balancing the condensate density and allowing us to accurately study the decay
of coherence across different coupling topologies. The result is a homogeneous
macroscopic, high energy, condensate lattice with coherence length exceeding
multiple lattice cells and near-diffraction limited emission. Generating such
large condensate lattices by using multiple excitation sources allows us to go
beyond the standard single-excitation source limit Wertz  _et al._ (2010);
Caputo  _et al._ (2018) overcoming beam size limitations, beam profile
inhomogeneities, and condensate fragmentation Kasprzak  _et al._ (2006). We
observe that the connectivity of the lattice significantly enhances the
system’s coherence length due to increased coherent coupling between adjacent
condensates.

##  II Results

We optically generate lattices of coupled polariton condensates using non-
resonant, pulsed and tightly-focused laser excitation spots for each
condensate node (see Methods for details). Polaritons generated at each
condensate node convert inherited potential energy into kinetic energy
resulting in a radially (ballistically) expanding polariton fluid from an
antenna like source Richard  _et al._ (2005); Töpfer  _et al._ (2020). When
two or more expanding condensates are brought together, interference effects
are revealed implying phase synchronisation between the condensation centres
and the emergence of a macroscopic order parameter Tosi  _et al._ (2012);
Cristofolini  _et al._ (2013); Ohadi  _et al._ (2016); Berloff  _et al._
(2017); Töpfer  _et al._ (2020); Alyatkin  _et al._ (2020).

The spatial coherence in a network (or lattice) of $N$ polariton condensates
generated under pulsed excitation is described by the integrated complex
coherence factor (see Appendix A)

$\tilde{\mu}_{ij}=\frac{\int\Gamma_{ij}(t)\mathrm{d}t}{\sqrt{\int\Gamma_{ii}(t)\mathrm{d}t\int\Gamma_{jj}(t)\mathrm{d}t}},\quad
i,j=1,...,N.$ (1)

Here, the correlation function
$\Gamma_{ij}(t)=\left\langle\psi_{i}(t)^{*}\psi_{j}(t)\right\rangle$ denotes
the mutual intensity of each pair of condensates averaged over many
realisations (pulses) of the system, and $\psi_{i}(t)$ is the complex-valued
amplitude of the $i\mathrm{th}$ condensate. While the modulus of the complex
coherence factor $|\tilde{\mu}_{ij}|\leq 1$ is a normalised measure for the
coherence between two condensate nodes its argument
$\tilde{\theta}_{ij}=\arg(\tilde{\mu}_{ij})$ represents their average phase-
difference.

###  A Condensate density stabilisation

We utilise a reflective liquid-crystal phase-only spatial light modulator
(SLM) to modulate the Gaussian excitation pump beam in the Fourier plane of
the optical excitation system to generate desired excitation pump spot
geometries at the focal plane of the microscope objective lens (see schematic
in Fig. 1(a)).

Figure 1: Closed-loop density stabilisation in polariton lattices. (a)
Schematic of the feedback loop that iteratively analyses the emission and
adjusts the pump profile to equalise the emission intensities of all
condensate nodes. (b) Integrated emission of 61 condensate nodes at
condensation threshold ($P=P_{\mathrm{thr}}$) in a triangular lattice
configuration (lattice constant $a=14.9\;\mathrm{\upmu m}$) without (blue
squares) and with (red circles) density stabilisation. Inset depicts the
condensate number indexing. (c,d) Recorded real-space photoluminescence and
(e,f) measured complex coherence factor $\tilde{\mu}_{1j}$ between the central
condensate node $1$ and each other condensate node $j$ without and with
density stabilisation, respectively. False colourscale and pseudo-spins (black
arrows) depict magnitude $|\tilde{\mu}_{1j}|$ and phase $\tilde{\theta}_{1j}$.
Data shown in (c-f) are extracted at a total pump power
$P=1.2P_{\mathrm{thr}}$. Scale bars in (c-f) correspond to $20\;\mathrm{\upmu
m}$. Abbreviations in (a): SLM - spatial light modulator, BS - beam splitter,
MO - microscope objective lens, MC - microcavity, LP - longpass filter, CCD -
Charge-coupled device.

Phase holograms are calculated using a modified version of the well-known
Gerchberg-Saxton (GS) algorithm Gerchberg and Saxton (1972), in which the
measured polariton photoluminescence (PL) of the condensate lattice feeds back
into the closed-loop sequence (see Appendix B for a full description of the
method). Iterative nonlinear adjustment of the excitation pump profile allows
us to create macroscopic lattices of $>100$ condensates with homogeneous
condensate node density ($\leq 1\%$ relative standard deviation (RSD)) across
the whole lattice.

In Fig. 1(b) we compare the measured emission of all condensate nodes in a
triangular lattice of 61 elements and lattice constant $a=14.9\;\mathrm{\upmu
m}$ when using $100$ iterations of the conventional GS algorithm (blue
squares) and the presented condensate density stabilisation method (red
circles), respectively. Unavoidable effects such as sample disorder, finite
accuracy of the GS algorithm, optical aberrations and missing translational
invariance due the finite size of the lattice Kalinin and Berloff (2018) all
contribute to a broad distribution of effective gain for each condensate node
and, thus, to a large spread $\approx 37\%$ (RSD) in the distribution of
condensate emission powers for the case of no condensate density
stabilisation. Active stabilisation of the lattice using the described closed-
loop sequence allows us to compensate for these detrimental elements and yield
a reduced spread $\approx 1\%$ (RSD), which is limited by experimental noise
in the system.

The near-field (or real space) polariton photoluminescence (PL) measured at
pump power $P=1.2P_{\mathrm{thr}}$, where $P_{\mathrm{thr}}$ is the system’s
condensation threshold pump level, is shown in Figs. 1(c) and (d) for both
excitation schemes. The presence of interference fringes with large visibility
in-between ballistically expanding condensates demonstrates mutual coherence
between nearest neighbour condensates. However, the lack of homogeneity in the
spatial distribution of interference fringes in Fig. 1(c) for the case of no
density stabilisation indicates a broad distribution of relative phase-
difference $\theta_{ij}$ between nodes. This is further confirmed by detailed
measurements of the integrated complex coherence factor $\tilde{\mu}_{ij}$
(see Appendix C and D for methods). In Figs. 1(e) and (f) we plot magnitude
$|\tilde{\mu}_{1j}|$ and phase $\tilde{\theta}_{1j}$ between the central
condensate node $1$ and each other condensate node $j$ using false-colour and
pseudo-spins (black arrows). We find an enhanced and isotropic spatial decay
of coherence $|\tilde{\mu}_{1j}|$ and larger homogeneity in relative phase-
differences $\tilde{\theta}_{1j}$ for the density stabilised polariton
lattice. The noticeable increase of phase-differences $\tilde{\theta}_{1j}$ \-
or analogous rotation of pseudo-spins - towards the edges of the lattice is an
expected finite-size effect due to a flux of particles escaping the system
Kalinin and Berloff (2018). In the ideal scenario of an infinite triangular
lattice, with homogeneous condensate occupation numbers, one retrieves a
homogeneous distribution of phase differences due to the system’s
translational invariance (see Appendix E).

Figure 2: Emission characteristics of a triangular polariton lattice in
reciprocal space. First Brillouin zones of the far-field emission for a
triangular polariton lattice of 61 condensates excited (a) without feedback,
(b) with feedback and for (c) a calculated ideal system of 61 coherent phase-
synchronised point sources. (d) Extracted profiles of the central Bragg peak
along $k_{x}=0$ for the far-field emission patterns shown in (a-c). (e)
Extracted FWHM of the central Bragg-peak as a function of total excitation
pump power. Data shown in (a,b,d) are recorded at pump power
$P=1.2P_{\mathrm{thr}}$.

The polariton condensate lattice is further probed by far-field (or reciprocal
space) measurements in analogy to time-of-flight experiments in cold atom
systems Hadzibabic  _et al._ (2006) (see Appendix C for methods). A zoom into
the first Brillouin zones of the recorded far-field emission pattern is shown
in Figs. 2(a) and (b) for the two lattice realisations with and without node
density stabilisation pumped at $P=1.2P_{\mathrm{thr}}$. For comparison, we
illustrate in Fig. 2(c) the calculated far-field emission of 61 superimposed
fully-coherent point-sources (wavefunctions) in triangular arrangement which
shows good agreement with the emission pattern of the stabilised condensate
lattice. We compare the width of the central far-field emission peak at $k=0$
for the three cases shown in Figs. 2(a-c) by plotting the corresponding
extracted intensity profiles along $k_{x}=0$ in Fig. 2(d). The lobe of the
calculated fully-coherent system (black line) represents the diffraction
limited interference peak. We extract the central lobe’s full-width-at-half-
maximum (FWHM) for varying pump power $P$ as shown in Fig. 2(e) and find near-
diffraction limited far-field emission for the density stabilised lattice (red
dash-dotted line), i.e. a minimum width at $P=1.17P_{\mathrm{thr}}$ which is
only $\approx 13\%$ larger than the diffraction limit of a fully-coherent
system (black horizontal line). Without active condensate density
stabilisation (blue dashed line) the peak width increases to $\approx 47\%$ at
the same pump power. The observed broadening of the far-field emission peak in
the physical system is a result of both reduced coherence $|\tilde{\mu}_{ij}|$
and non-homogeneous phase distribution $\tilde{\mu}_{ij}$. The system’s pump-
power dependencies are summarised in Appendix F.

###  B Coherence vs. dimensionality

In this section we begin by investigating the coherence between two ballistic
polariton condensates. Previous studies of this system have revealed
periodically alternating synchronisation patterns of in-phase
($\tilde{\theta}_{12}=0$) and anti-phase ($\tilde{\theta}_{12}=\pi$) states
with increasing condensate separation distance $d_{12}$ Ohadi  _et al._
(2016); Töpfer  _et al._ (2020). While there exist separation distances at
which the emission of the coupled condensate system is not single-mode but
exhibits two or more modes of both even and odd parity states Ohadi  _et al._
(2016); Töpfer  _et al._ (2020), in the following, we focus on single-mode
realisations only. Under this condition the integrated complex coherence
factor $\tilde{\mu}_{12}$ is a measure for the system’s average coherence
properties. In Appendix G we detail on the time-resolved coherence build-up
$\mu_{12}(t)$ of a polariton dyad.

Figure 3: Ballistic expansion, coupling and interference of two polariton
condensates. Recorded (a,d) real-space and (b,c) far-field PL of two
condensates with separation distances $d_{12}=12.7\;\mathrm{\upmu m}$ and
$d_{12}=89.3\;\mathrm{\upmu m}$, respectively. Corresponding far-field
interference patterns after masking of the emission in real space to block all
emission outside the $2\;\mathrm{\upmu m}$ FWHM of each condensate nodes are
shown in (c) and (f). (g) Distance dependence of the extracted integrated
complex coherence factor $|\tilde{\mu}_{12}|$, while keeping the excitation
pump power constant at $P=1.2P_{\mathrm{thr}}$, where $P_{\mathrm{thr}}$ is
the measured threshold pump power at a separation distance of
$d_{12}=12.7\;\mathrm{\upmu m}$. Red curve is a Gaussian fit (Eq. (2)) to the
experimental data points. Inset shows the mutual coherence obtained by
numerical simulations. False colour scale in (f) applies to (b,c,e,f) in
linear scale and to (a,d) in logarithmic scale saturated below $10^{-4}$ of
the maximum count rate. Scale bars in (a,d) and (b,c,e,f) correspond to
$10\;\mathrm{\upmu m}$ and $1\;\mathrm{\upmu m}^{-1}$, respectively.

Recorded near-field and far-field PL of two condensates separated by
$d_{12}=12.7\;\mathrm{\upmu m}$ and pumped equally at $P=1.2P_{\mathrm{thr}}$
are shown in Figs. 3(a) and (b). The interference patterns in both emission
images reveal anti-phase synchronisation between the two condensates
($\tilde{\theta}_{12}=\pi$) due to the even number of interference fringes.
Far-field emission consists predominantly of PL at large in-plane wavevector
$k=1.8\;\mathrm{\upmu m}^{-1}$ demonstrating the ballistic expansion of both
condensates due to the strong repulsive interactions between polaritons and
the pump-induced exciton reservoirs localised at each condensate center
Wouters  _et al._ (2008). In Fig. 3(c) we show the recorded far-field
interference pattern when spatially filtering the central $2\;\mathrm{\upmu
m}$ FWHM of each condensate. Analysis of this double hole interference pattern
reveals a large coherence factor $|\tilde{\mu}_{12}|=0.9$. Increasing the
separation distance $d_{12}$ while keeping the the pump power $P$ constant
leads to a decay of coherence $|\tilde{\mu}_{12}|$ between the two condensate
centres such that at $d_{12}=89.3\;\mathrm{\upmu m}$ no interference pattern
can be found in near-field an far-field emission as shown in Figs. 3(d-f).

The dependency of the measured coherence $|\tilde{\mu}_{12}|$ on the pump spot
separation distance $d_{12}$ is shown in Fig. 3(g) as blue circles and is
fitted with a Gaussian decay

$|\tilde{\mu}_{12}(d_{12})|=\exp{\left\\{-\frac{\pi}{4}\left(\frac{d_{12}}{L_{c}}\right)^{2}\right\\}}.$
(2)

The fit parameter $L_{c}=24.9\pm 0.4\;\mathrm{\upmu m}$ represents a measure
for the length over which synchronisation of the two-condensate system is
possible, and we denote it as the effective coherence length,

$L_{c}=\int_{0}^{\infty}|\tilde{\mu}_{12}(x)|\mathrm{d}x.$

We reproduce the experimentally measured decay of mutual coherence between two
coupled polariton condensates through 2D numerical simulation using the
generalised stochastic Gross-Pitaevskii equation (GPE) shown in Fig. 3(g) as
yellow squares (see Methods). The dependency of coherence $|\tilde{\mu}_{12}|$
as a function of pump power $P$ between two coupled condensates with
separation distance $d_{12}=12.7\;\mathrm{\upmu m}$ is depicted in the inset
of Fig. 3(g) for both experiment (blue circles) and numerical simulation
(yellow squares). An increase of coherence between the two condensate nodes
arising at condensation threshold is followed by a drop of
$|\tilde{\mu}_{12}|$ for larger pump power $P>1.2P_{\mathrm{thr}}$. This
decrease in $|\tilde{\mu}_{12}|$ is largely attributed to the transition of
the system into multi-mode operation Ohadi  _et al._ (2016) yielding reduced
visibility in time-integrated measurements. We point out that the system
realisations which are shown in the distance-dependence in Fig. 3(g) exhibit
single-mode emission only. The decrease in coherence with increasing
separation distance is caused by the spatial decay of the wavefunction of
ballistically expanding polariton condensates Töpfer  _et al._ (2020), i.e. a
reduced coupling strength between the two condensate nodes. In Appendix H we
compare our results to a dyad under continuous wave excitation.

Figure 4: Finite 1D lattice of ballistically coupled polariton condensates.
Recorded (a) near-field and (b) far-field PL for 11 condensates arranged in a
chain with lattice constant $a=12.1\;\mathrm{\upmu m}$. (c) Correlation matrix
illustrating measured mutual coherence factor $|\tilde{\mu}_{ij}|$ and phase
factor $\tilde{\theta}_{ij}=\arg{(\tilde{\mu}_{\ij})}$ for each pair of
condensates $i\neq j$. Condensate indexing is shown in (a). (d) Spatial decay
of the averaged mutual coherence factors $|\tilde{\mu}_{ij}|$ with GPE
simulation shown in inset. Error bars represent the standard deviation of the
set of coherence factor values for each separation distance $d_{ij}$. Red
curve is an exponential fit to the experimental data points. The pump power is
kept at $P=1.2P_{\mathrm{thr}}$ for all data shown in (a-d). Scale bars in (a)
and (b) correspond to $20\;\mathrm{\upmu m}$ and $1\;\mathrm{\upmu m}^{-1}$,
respectively.

Next, we increase the number of condensates by investigating a linear chain of
11 equally spaced condensates with nearest-neighbours (NN) distance
$a=12.1\;\mathrm{\upmu m}$. The resultant interference patterns in both real
and reciprocal space shown in Figs. 4(a) and (b) indicate anti-phase
synchronisation, i.e. phase differences between any pair $\\{i,j\\}$ of
condensates according to

$\tilde{\theta}_{ij}=\pi\cdot(i-j)\mod\pi,$ (3)

with condensate indexing shown in Fig. 4(a). Such order can be said to be
“antiferromagnetic”.

We measure and analyse the far-field interference between each pair of
condensates $i\neq j$. Extracted magnitude $|\tilde{\mu}_{ij}|$ and phase
$\tilde{\theta}_{ij}$ of the integrated complex coherence factor are
illustrated in matrix form in Fig. 4(c) with row and column indices $i,j$
denoting the pair of condensates, where we make use of the hermiticity of the
correlation matrix $\tilde{\mu}_{ij}=\tilde{\mu}^{*}_{ji}$. In agreement with
Eq. (3) we confirm anti-phase synchronisation between NNs (antiparallel
pseudo-spins). For large condensate pair distances, $|i-j|\gg 1$, the
coherence factor $\tilde{\mu}_{ij}$ decays in magnitude and its phase deviates
from Eq. (3) indicating loss of long-range antiferromagnetic order. In Fig.
4(d) we show the decay of coherence $|\tilde{\mu}_{ij}|$ in the chain as a
function of condensate spacing $d_{ij}$ (i.e., lattice constant $a$ is fixed).
The data are fitted with a single exponential decay (red curve)

$|\tilde{\mu}_{ij}|=\exp{\left\\{-d_{ij}/L_{c}\right\\}},$ (4)

which yields a coherence length $L_{c}=35\pm 1\;\mathrm{\upmu m}$.

In a next step, we increase the system size to a 2D square lattice comprising
121 condensates with lattice constant $a=12.1\;\mathrm{\upmu m}$. The
interference patterns in near-field and far-field PL which are shown in Figs.
5(a) and (b) reveal anti-phase synchronisation in analogy to the 1D system in
Fig. 4. We measure $\tilde{\mu}_{1j}$ between the central condensate (index
$1$) and all other condensate nodes $j=2,...,121$ plotted in Fig. 5(c). The
observed antiferromagnetic order ($\tilde{\mu}_{1j}\approx\pm 1)$ reduces with
distance from the centre seen from the pseudo-spin rotation towards the edges
in Fig. 5(c) due to flux of particles escaping the lattice. However, coherence
$|\tilde{\mu}_{1j}|$ between the most central condensate and any other
condensate of the lattice does not drop below $0.3$ for separation distances
as large as $d_{1j}=86\;\mathrm{\upmu m}$ towards the corners of the lattice.
The spatial decay of coherence $|\tilde{\mu}_{1j}|$ is shown in Fig. 5(d) and
fitted with an exponential function (Eq.(4), red curve) with coherence length
$L_{c}=87\pm 1\;\mathrm{\upmu m}$.

Figure 5: Finite two-dimensional square lattice of ballistically coupled
polariton condensates. Recorded (a) near-field and (b) far-field PL for 121
condensates arranged in a square lattice with lattice constant
$a=12.1\;\mathrm{\upmu m}$. Spatial correlation map showing mutual coherence
factor $|\tilde{\mu}_{1j}|$ and phase factor
$\tilde{\theta}_{1j}=\arg{(\tilde{\mu}_{1j})}$ between the most central
condensate (index 1) and each other condensate (index $j$). (d) Spatial decay
of the averaged mutual coherence factors $|\tilde{\mu}_{1j}|$ with GPE
simulation shown in inset. Error bars represent the standard deviation of the
set of coherence factor values for each separation distance $d_{ij}$. Red
curve is an exponential fit to the experimental data points. The pump power is
kept at $P=1.2P_{\mathrm{thr}}$ for all data shown in (a-d). Scale bars in (a)
and (b) correspond to $20\;\mathrm{\upmu m}$ and $1\;\mathrm{\upmu m}^{-1}$,
respectively.

To compare the obtained results, we summarise in Fig. 6(a) the spatial decay
of coherence for four different types of networks; two coupled condensates, 1D
chain, and 2D square and triangular lattices. For the lattices the data points
for different separation distances $d_{ij}$ correspond to different pairs of
condensate nodes $\\{i,j\\}$. For the polariton dyad the physical separation
distance between the two condensate nodes is changed. In all cases the pump
power was chosen such as to maximise the coherence $|\tilde{\mu}_{12}|$
between a pair of NNs in the system. In Fig. 6(c) we illustrate the extracted
(effective) coherence length $L_{c}$ versus the average number of NNs in each
system. It is apparent that the coherence length $L_{c}$ is increasing for
systems with larger connectivity which we explain in the following.
Interestingly, the mutual coherence between NNs in the 1D chain ($d_{ij}=a$)
is noticeably lower than in the other systems. Perhaps most surprisingly, the
chain’s coherence is lower than the dyad for small $d_{ij}$ and then overtakes
it around $d\gtrsim 25\;\mathrm{\upmu m}$.

Figure 6: Comparison of coherence between ballistically coupled polariton
condensates in different networks presented in Figs. 1, 3, 4 and 5. (a)
Spatial decay of coherence for a two condensate system of varying size (grey
triangles), a 1D chain (blue squares), a square lattice (red squares) and a
triangular lattice (black diamonds). Curves represent Gaussian and exponential
fits to the experimental data points. (b) Excitation pulse energy dependence
of the mutual coherence factor $|\tilde{\mu}_{12}|$ between the central-most
condensate (index $1$) and one of its NNs (index $2$) for the 1D and 2D
periodic structures, and between two condensates at
$d_{12}=12.7\;\mathrm{\upmu m}$ for the case of a dyad. The pump powers used
in (a) are marked with vertical arrows and are given by $P\approx
1.2P_{\mathrm{thr}}$, where $P_{\mathrm{thr}}$ is the threshold pump power for
each system, respectively. (c,d) Coherence length $L_{c}$ and average
excitation pulse energy $E_{\mathrm{thr}}$ per condensate node at threshold
versus the average number of NNs $Z$. (e,f) Calculated Lyapunov exponents in
$\Gamma\to X$ crystal momentum space ($q$-space) around the $\Gamma$-point
($qa=0$) belonging to a square condensate lattice with small (high) particle
numbers given by yellow (black) curves.

We attribute the origin of increased spatial coherence for condensate networks
with larger connectivity to their reduced population of reservoir excitons as
we argue in the following. In Fig. 6(b) we compare the dependence of coherence
$|\tilde{\mu}_{12}|$ between a pair of NNs as a function of the average
excitation pulse energy per condensate in each system. We observe a reduction
of threshold pulse energy $E_{\mathrm{thr}}$ per condensate node for networks
with larger connectivity (larger number of NNs), which is shown in Fig. 6(d)
and attributed to the increased gain given by ballistic exchange of particles
between neighbouring nodes. A lower pump energy generates a smaller reservoir
of excitons for each condensate node and since reservoir-condensate
interactions play a dominant role in condensate dephasing Wouters and Savona
(2009), a lower threshold pump energy will generally result in increased
coherence properties.

By discretising the system into a set of interacting condensates Töpfer  _et
al._ (2020) one can show that threshold pump power and coherence length scale
approximately linear with the number of condensate NNs. Let us consider the
steady state (without phase frustration) of frequency $\nu$ such that the
discretised Gross-Pitaevskii equation can be written,

$\nu\psi_{n}=\left[V+\frac{i}{2}\left(\frac{P}{n_{s}+|\psi_{n}|^{2}}-\gamma\right)\right]\psi_{n}+(\epsilon+i\kappa)\sum_{\langle
nm\rangle}\psi_{m},$ (5)

Here, $\psi_{n}$ describes the phase and amplitude of the $n$th with sum over
NNs, the potential $V$ describes blueshift coming from the condensate
interactions with itself and its noncondensed particle background, $P$ is the
laser power, $n_{s}$ is the condensate saturation density, and $\gamma$ is the
polariton lifetime. The non-Hermitian interaction between condensates is
captured with $\epsilon$ and $\kappa$. The condensation threshold (the lasing
mode) belongs to the mode with the lowest particle decay and its background
reservoir of noncondensed particles responsible for dephasing is written,

$n_{x}^{\text{(thr)}}=n_{s}-\frac{Z\kappa}{R},$ (6)

where $Z$ is the number of NNs and $R$ is the scattering rate of reservoir
particles into the condensate. The density of the reservoir
$n_{x}^{\text{(thr)}}$ is proportional to the amplitude of random fluctuations
experienced by the condensate (see Methods) and therefore its coherence
properties will depend on the number of neighbours $Z$. The first order
spatial correlation function can be written,

$\mu_{nm}\propto\exp{[-f_{\mathbf{k},\omega}(\mathbf{r}_{nm})]},$ (7)

where $\mathbf{r}_{nm}$ is the distance vector between condensates $n$ and $m$
and $f_{\mathbf{k},\omega}$ is the phase-phase correlation function which
determines the condensate’s long range coherence properties Szymańska  _et
al._ (2006). It is beyond the scope of this study to calculate
$f_{\mathbf{k},\omega}$ given the complicated spatial density of the
condensate. At the length scale of the experiment we observe an exponential
decay of spatial correlations in all cases (except for the polariton dyad)
shown in Fig. 6(a) and therefore we will assume that
$f_{\mathbf{k},\omega}=|\mathbf{r}_{nm}|/L_{c}$ holds. Assuming frequency
independent noise processes the coherence length $L_{c}$ becomes inversely
proportional to the strength of the noise Szymańska  _et al._ (2006); Roumpos
_et al._ (2012) and one obtains,

$L_{c}\propto\frac{1}{\gamma+Rn_{s}-Z\kappa}\simeq\frac{1}{\gamma+Rn_{s}}\left[1+\frac{Z\kappa}{\gamma+Rn_{s}}\right].$
(8)

The expected linear dependencies of both threshold pump power
$P_{\mathrm{thr}}\propto n_{x}^{\text{(thr)}}$ (Eq. (6)) and coherence length
$L_{c}$ (Eq. (8)) are both indicated as gray lines in Figs. 6(c,d). We point
out that similar results (increasing coherence with larger connectivity and
dimensionality in networks of coupled elements) have been observed in other
technological platform such as in arrays of coupled VCSELs Golshani  _et al._
(1999); Pier  _et al._ (2000), micromechanical oscillator arrays Zhang  _et
al._ (2015) as well as for coupled fiber lasers Fridman  _et al._ (2010).

It is instructive to investigate the condensate’s dispersion of elementary
excitations (Lyapunov exponents) through the GPE which directly relates to the
behaviour of space-time correlations in the system. For simplicity we will
focus on the case of a 2D square lattice such as displayed in Fig. 5 in a
steady state (continuous wave excitation). For our calculation of the Lyapunov
exponents we will assume that we are working in the the bulk of the condensate
and therefore the system is taken to have discrete translational invariance.
This allows us to apply Bloch’s theorem to the standard Bogoliubov treatment
and solve the Lyapunov exponents in the reduced Brillouin zone of the
condensate-pump lattice (see Appendix I). In Fig. 6(e,f) we plot the Lyapunov
exponents around the $\Gamma$-point of the reduced Brillouin zone for two
different numbers of particles in the condensate. Yellow and black coloured
curves correspond to $N\approx 100$ and $\approx 800$ particles in the
condensate unit cell respectively for a lattice constant of $a=12$ $\upmu$m.
$\text{Im}{(\lambda)}$ corresponds to oscillatory evolution of the
fluctuations, showing the purely diffusive Goldstone branch
[$\text{Im}{(\lambda)}=0$] at lower momenta. $\text{Re}{(\lambda)}$
corresponds to the decay (growth) rate of fluctuations when negative
(positive) valued. In the top branch of the pitchfork, we observe that
$\text{Re}{(\lambda)}$ becomes more negative for a given wavevector as more
particles are in the condensate (black curve). This is in agreement with our
experiment where coherence increases with pump power (particle number grows
and decay of fluctuations increases) up to the point where single-mode
(stationary) behaviour is lost. Our results are similar to those obtained for
spatially uniform systems Szymańska  _et al._ (2006); Wouters and Carusotto
(2007); Chiocchetta and Carusotto (2013) underlining that spatial details of
the condensate structure are not relevant to long-wavelength fluctuations.

##  III Discussion

We have demonstrated that the coherence properties in lattices of polariton
condensates are enhanced by balancing the condensate emission across the
system using a closed-loop feedback scheme to adjust the excitation pump
geometry. This scheme reduces the effects of optical aberrations in the
experimental system, as well as counteracts mode localisation due to sample
non-uniformities Pier and Kapon (1997) and non-homogeneous gain distribution
across the coupled condensate network. While actively controlling the
condensate lattice uniformity we have accurately determined phase and
coherence between any pairs of condensates in different types of networks
ranging from two condensate-systems to 1D and 2D periodic structures. The
dynamics and synchronisation of coupled non-linear elements is critically
influenced by the underlying coupling topology Watts and Strogatz (1998).
Here, we have shown that an increase in connectivity, i.e. number of NNs,
significantly reduces the operational pump power per node and increases the
coherence length demonstrating a promising route to polaritonic devices with
networks of many coupled condensates with potential application in simulators
and optical based computation. The presented measurements and techniques
provide a deeper understanding on the coherence properties of coupled light-
matter wave fluids in low dimensional quantum systems, and qualify for other
open (dissipative) networks such as laser arrays and photon condensate
lattices. Furthermore, our analysis of the lattice condensate fluctuations
reveals similar long-wavelength dispersions to those of uniform systems which
are regarded as the ideal case from a theoretical point of view. This suggests
access to fundamental long-wavelength physics belonging to uniform systems by
designing instead an extended condensate lattice.

##  IV Acknowledgements

The authors acknowledge the support of the Skoltech NGP Program (Skoltech-MIT
joint project), the UK’s Engineering and Physical Sciences Research Council
(grant EP/M025330/1 on Hybrid Polaritonics), and the RFBR projects No.
20-52-12026 (jointly with DFG) and No. 20-02-00919. Y.G.R. acknowledges
support by CONACYT (Mexico) Grant No. 251808 and by PAPIIT-UNAM Grant No.
IN106320.

##  V Methods

###  A Experimental methods

We utilise a strain-compensated $2\lambda$ GaAs microcavity with embedded
InGaAs quantum wells and a quality factor $Q\approx 12\;000$ Cilibrizzi  _et
al._ (2014) . We continuously cool the microcavity sample in a cold-finger
cryostat ($T\approx 6\;\mathrm{K}$) and operate at a negative cavity-detuning
of $\Delta\approx-5\;\mathrm{meV}$ resulting in a lower polariton mode at
$\lambda\approx 858\;\mathrm{nm}$ for zero in-plane momentum $|\mathbf{k}|=0$.
We use a circularly polarised non-resonant blue-detuned ($\lambda\approx
800\;\mathrm{nm}$) pulsed laser (pulse duration $\approx 150\;\mathrm{fs}$,
repetition rate $80\;\mathrm{MHz}$). For each condensate node the laser is
focused onto the microcavity sample to a beam waist of approximately
$2\;\mathrm{\upmu m}$ via a $0.4\;\mathrm{NA}$ microscope objective. Polariton
photoluminescence (PL) and the excitation laser are measured by imaging the
cavity emission or reflected light of the microcavity sample onto a CCD sensor
with integration times in the range of $10$-$100\;\mathrm{ms}$, which
corresponds to averaging measurements over approximately one million
realisations of the system. By using a longpass or shortpass filter in front
of the camera we can selectively choose to measure either the spatial geometry
of the polariton emission or the excitation laser.

Experimental details of the density stabilisation technique for polariton
condensate lattices are given in Appendix B. Methods for the measurement of
the mutual complex coherence factor $\mu_{ij}$ between any pair of condensate
nodes are described in Appendices C and D.

###  B Error analysis

Error bars for the absolute values of the complex coherence factor
$|\mu_{ij}|$ represent the confidence interval of the fitted parameter and are
extracted as described in Ref. Singer  _et al._ (2012). We fit the double hole
interference pattern (Eq. 13) to the experimental interference pattern
$I_{1+2}$ with three fit parameters $|\mu_{ij}|$, $|\mathbf{d}_{12}|$ and
$\theta_{12}$. For each measurement we define the quality of the fitted
interference curve $I_{\mathrm{fit}}$ by the dimensionless quantity

$R=\frac{\sum_{k_{||}}\left[I_{\mathrm{fit}}(k_{||})-I_{1+2}(k_{||})\right]^{2}}{\sum_{k_{||}}\left[I_{1+2}(k_{||})\right]^{2}}$

where summation is performed over all sample points $k_{||}$. Next, we keep
all fit parameters constant except the coherence factor $|\mu_{ij}|$ and
recalculate the quality factor $R$ for varying $|\mu_{ij}|$. The confidence
interval of $|\mu_{ij}|$ is determined as the range in which the quality
factor $R$ increases by $100\%$ from its minimum value.

###  C Numerical Simulations

A stochastic dissipative Gross-Pitaevskii equation describes the polariton
condensate order parameter $\Psi(\mathbf{r},t)$ coupled to a rate equation
describing the density of a background excitonic reservoir $n_{x}$ feeding
particles into the condensate. The reservoir itself is sustained by a decaying
population of excited electron-hole pairs
$n_{c}(\mathbf{r},t)=P(\mathbf{r})e^{-(W+\Gamma_{c})t}$ which we assume are
generated instantaneously by the sub-picosecond nonresonant Gaussian shaped
pump $P(\mathbf{r})$,

$\displaystyle i\frac{\partial\Psi}{\partial t}$
$\displaystyle=\frac{1}{2}\left[-\frac{\hbar\nabla^{2}}{m}+V+\alpha|\Psi|^{2}+i\left(Rn_{x}-\gamma\right)\right]\Psi+\frac{dW}{dt},$
(9a) $\displaystyle\frac{\partial n_{x}}{\partial t}$
$\displaystyle=-\left(\Gamma_{R}+R|\Psi|^{2}\right)n_{x}+Wn_{c}.$ (9b)

Here, $m$ is the effective mass of a polariton in the lower dispersion branch,
$\alpha$ is the interaction strength of two polaritons in the condensate, $g$
is the polariton-reservoir interaction strength, $R$ is the rate of stimulated
scattering of polaritons into the condensate from the active reservoir,
$\gamma$ is the polariton condensate decay rate, $\Gamma_{R,c}$ are the decay
rates of the reservoir excitons and electron-hole pairs respectively, $W$ is
the conversion rate between the reservoir excitons and electron hole-pairs,
and $V(\mathbf{r},t)=g(n_{x}(\mathbf{r},t)+n_{c}(\mathbf{r},t))$ is the pump-
induced potential with an effective reservoir-condensate interaction strength
$g$. Here we have introduced a Gaussian white noise term $dW$ based on a Monte
Carlo technique in the truncated Wigner representation Wouters and Savona
(2009). The correlators are written $\langle
dW(\mathbf{r}_{i})dW(\mathbf{r}_{j})\rangle=0$ and $\langle
dW^{*}(\mathbf{r}_{i})dW(\mathbf{r}_{j})\rangle=(\Gamma+Rn_{x})\delta(\mathbf{r}_{i}-\mathbf{r}_{j})dt/2\Delta
A$ where $\Delta A$ is the cell area of the spatial grid. The coherence factor
(Eq. II) is then calculated by averaging over multiple realisations of
condensate formation. Good agreement between simulation and experimental data
is obtained by only adjusting the values of $g$ and
$\text{max}[P(\mathbf{r})]$ between different systems. The reason $g$ is taken
as an tuneable parameter is due to the fact that $V$ can possess a more
complicated sublinear dependence in Eq. (9a) as pointed out previously Gao
_et al._ (2015). Parameters were set to: $\Gamma=1/5.5$ ps-1,
$\Gamma_{R}=2\Gamma$, $\Gamma_{c}=0.0055\Gamma$, $W=0.275\Gamma$,
$\alpha=0.01$ ps-1 $\upmu$m2, $R=0.05\alpha$, $m=0.32$ meV ps2 $\upmu$m-2, and
$g=\\{0.0008,\ 0.0011,\ 0.0022,\ 0.0016\\}\alpha$ for dyad, chain, square, and
triangle respectively.

##  APPENDIX A COHERENCE IN POLARITON LATTICES

In our study of synchronisation in networks of coupled polariton condensates
we are interested in the correlations between pairs of condensates which are
described by their mutual coherence function

$\Gamma_{ij}(t)=\left\langle\psi_{i}(t)^{*}\psi_{j}(t)\right\rangle,$

where the brackets denote ensemble averaging and $\psi_{i}(t)$ is the complex-
valued amplitude of the $i$th condensate node. The expectation value of the
condensate’s particle number (occupation) is given by the diagonal elements
$\Gamma_{ii}(t)$, and the presence of non-zero off-diagonal elements
$|\Gamma_{ij}(t)|>0$ for $i\neq j$ indicates long-range order in the coupled
condensate network. For nonstationary systems (such as for polariton
condensates generated under pulsed excitation) the coherence function
$\Gamma_{ij}(t)$ depends on time $t$. We note that pulsed excitation ($\approx
150\;\mathrm{fs}$ pulse width) of ballistically coupled condensate nodes leads
to polariton photoluminescence with typical emission signal full-width-at-
half-maximum (FWHM) of $18\;\mathrm{ps}$ (see Appendix G). A normalised form
of the mutual coherence function is given by the complex coherence factor

$\mu_{ij}(t)=\frac{\Gamma_{ij}(t)}{\sqrt{\Gamma_{ii}(t)\Gamma_{jj}(t)}},$

where for a fully coherent pair of condensates $|\mu_{ij}|$ = 1. Considering
measurements with time integration over many pulses we define the integrated
complex coherence factor

$\tilde{\mu}_{ij}=\frac{\int\Gamma_{ij}(t)\mathrm{d}t}{\sqrt{\int\Gamma_{ii}(t)\mathrm{d}t\int\Gamma_{jj}(t)\mathrm{d}t}},$

as an averaged measure for the mutual coherence properties of condensates in
polariton networks. Throughout this manuscript time-integrated variables are
marked with the tilde diacritic. For stationary ergodic systems (such as
polariton condensates under constant pumping) the complex coherence factor
$\mu_{ij}(t)$ does not depend on time and time-integrated measurements of
correlations fully determine the coherence properties, i.e.
$\mu_{ij}(t)=\tilde{\mu}_{ij}$. In case of pulsed excitation of polariton
condensates the modulus of the time-integrated factor $|\tilde{\mu}_{ij}|$
represents a lower bound for the maximum value of the complex coherence factor
$|\mu_{ij}(t)|$. We note that the argument of the integrated complex coherence
factor represents the average condensate node phase difference
$\tilde{\theta}_{ij}=\arg(\tilde{\mu}_{ij})$ and the integrated particle
number $\int\Gamma_{ii}(t)\mathrm{d}t$ of each condensate node is proportional
to its measurable average emission intensity $\tilde{I}_{i}$.

##  APPENDIX B FEEDBACK SCHEME

In this manuscript we present networks of ballistically coupled polariton
condensates for which the pumping geometry has been adjusted using a
reflective SLM (see Fig. 1) such as to equalise the emission power of each
condensate node at condensation threshold. The algorithm used to calculate the
holograms which are applied to the SLM is a modified version of the Gerchberg-
Saxton (GS) algorithm Gerchberg and Saxton (1972) and takes into account
feedback from the recorded PL. This iterative algorithm is similar to the
schemes presented in Refs. [Nogrette  _et al._ , 2014,Tamura  _et al._ , 2016]
applied to trapping of cold atoms. In our case the recorded PL of the
condensates is not directly proportional to the pump power of each node but
depends on the coupling-topology of the network, i.e. ballistic coupling
between spatially separated condensate nodes affects the gain of each
condensate.

The algorithm is schematically illustrated in Fig. 7(a) and begins with the
initialisation of a 2D complex field of amplitude $A^{(0)}(\mathbf{r})$ and
phase $\phi^{(0)}(x,y)$ representing the amplitude of the pump laser and the
phase hologram in pixel coordinates of the SLM, respectively. As an initial
guess we assume a constant phase pattern $\phi^{(0)}=0$ and a Gaussian field
amplitude, whose width matches the pump laser width. In each iteration step
$n$ propagation of the complex field $A^{(0)}e^{i\phi^{(n)}(\mathbf{r})}$ from
the SLM plane to the focal plane is computed by a fast Fourier transform (FFT)
yielding a complex-valued target field with amplitude
$\tilde{A}^{(n)}(\mathbf{r})$ and phase
$\phi_{\mathrm{GS}}^{(n)}(\mathbf{r})$. We add the output phase pattern
$\phi_{\mathrm{GS}}^{(n)}$ with a device dependent correction term $\phi_{c}$,
which can also take into account optical aberrations in the experimental
setup, and a blazed grating $\phi_{g}$ to offset the modulated laser beam from
its unmodulated reflected part. The sum of these terms modulo the $2\pi$
bitlevel is applied to the SLM and defines the pump laser geometry. In a next
step we record the near-field photoluminescence of the polariton system
excited by the current pump laser geometry and extract the emission power
$I_{\mathrm{C}}^{(n)}$ of each condensate by integrating the recorded signal
within the FWHM ($\approx 2\;\mathrm{\upmu m}$) of each node. We update the
target power $I_{\mathrm{T}}^{\mathrm{(n+1)}}$ for the next iteration step of
the algorithm for each condensate node utilising the non-linear function

$I^{(n+1)}_{T}=\frac{I^{(n)}_{T}}{1+\epsilon\left(I^{(n)}_{C}/\langle
I^{(n)}_{C}\rangle-1\right)},$ (10)

where $\epsilon$ is an adjustable feedback parameter and $\langle
I^{(n)}_{C}\rangle$ denotes the mean of the measured condensate node emission
distribution. The non-linear mapping given by Eq. (10) lowers (raises) the
target pump power $I^{(n+1)}_{T}$ of nodes with measured emission power
$I^{(n)}_{C}$ larger (lower) than the mean value $\langle I^{(n)}_{C}\rangle$.
The parameter $\epsilon$ controls the speed of the feedback loop but cannot be
set too large to avoid destabilisation of the algorithm. For $\epsilon=0$ the
algorithm represents the conventional form of the GS algorithm with no
iterative adjustment of target spot intensities, i.e.
$I^{(n+1)}_{T}=I^{(n)}_{T}$. Next, a complex-valued target field with
amplitude $\sqrt{I_{\mathrm{T}}}(\mathbf{r})$ and phase
$\phi_{\mathrm{GS}}^{(n)}(\mathbf{r})$ is constructed. Hereby, the field
amplitude $\sqrt{I_{\mathrm{T}}}(\mathbf{r})$ consists of the superposition of
Gaussian spots representing the position of each condensate node. The
amplitude of each Gaussian peak is given by the square root of the updated
spot power, i.e. $\sqrt{I_{\mathrm{T}}^{\mathrm{(n+1)}}}$ and the width is
related to the width of the Gaussian profile of the pump laser by means of a
Fourier transform. Back-propagation of this complex field to the SLM plane is
computed by an inverse fast Fourier transform (IFFT) yielding a phase pattern
$\phi^{(n+1)}$ which replaces the initial phase pattern for the next iteration
step, i.e. $\phi^{(n+1)}\rightarrow\phi^{(n)}$. The amplitude field at the SLM
plane is kept constant as the initial field $A^{(0)}$ for each iteration step.
We quit execution of the algorithm once the measured spot power distribution
$I_{\mathrm{C}}^{(n)}$ shows a spread which is smaller than a certain
threshold level (usually $1\%\;\mathrm{(RSD)}$ reached within
$<100\;$iterations).

Figure 7: (a) Schematic of the modified iterative Gerchberg-Saxton algorithm
including feedback of the polariton photoluminescence (PL) to stabilise the
node density in networks of coupled polariton condensates. (b,c) Recorded
laser pump profile and (d,e) real space PL at condensation threshold without
and with feedback, respectively.

We find good results for stabilising the intensity distribution in
ballistically coupled polariton networks ($\mathrm{RMS}\leq 1\%$) pumped at
condensation threshold $P\gtrsim P_{\mathrm{thr}}$ by choosing $\epsilon$ to
be in the range $10^{-2}-10^{-1}$. Stabilisation at larger pump powers $P\gg
P_{\mathrm{thr}}$ is impeded using the presented technique - which operates at
a typical rate of $60\;\mathrm{Hz}$ \- because of stronger nonlinear effects
resulting in unstable regimes Wouters and Carusotto (2007); Bobrovska  _et
al._ (2014); Baboux  _et al._ (2018) and/or excitement of polaritons to higher
energy states Wertz  _et al._ (2010), with subsequent nonstationary condensate
dynamics on a picosecond timescale.

In Fig. 7(b) we show the recorded laser (pump) profile for a triangular
lattice of 61 nodes (lattice constant $a=14.9\;\mathrm{\upmu m}$) using a
phase hologram computed with the conventional GS algorithm
($\epsilon=0;100\;\mathrm{iterations}$). Although the target image is a laser
pattern with homogeneous spot intensities, due to limited accuracy of the
algorithm as well as unavoidable optical aberrations due to misalignment and
device imperfections, the resultant experimentally obtained pump spot
intensity distribution deviates from the homogeneous target. We measure a
spread of $\approx 17\%$ (RSD) for the distribution of pump spot powers. Using
the same pump profile to excite a triangular lattice of polariton condensates
at condensation threshold ($P\approx P_{\mathrm{thr}}$) we record the near-
field PL as shown in Fig. 7(d) and find an increased spread of condensate node
densities of $\approx 37\%\;\mathrm{(RSD)}$. While sample disorder can have an
impact on the spatially dependent condensation threshold pump power and
emission intensity for each condensate node, the main reason for the increased
spread in the condensate emission power distribution is due to the system’s
nonlinear input-output-power characteristics and the finite size of the
system. Nodes that are positioned at the edge of the condensate network are
coupled to fewer nearest-neighbours (NN) than nodes in the bulk. Since
coherent coupling between ballistically expanding polariton condensates
results in a reduced condensation threshold Cristofolini  _et al._ (2013),
nodes coupled to fewer NNs will typically have a higher threshold pump power
than nodes with a greater number of NNs (see section B). As a result,
condensate nodes in the bulk of the network generally have a larger occupation
number (emission) than condensate nodes at the edges when pumped with the same
excitation pump power Kalinin and Berloff (2018). This is a consequence of
polariton waves radiating to the outside continuum and escaping the network
which can be regarded as an effective flux-induced potential Ostrovskaya  _et
al._ (2013).
In Figs. 7(c,e) we show the recorded pump laser and corresponding near-field
polariton emission profiles at threshold $P\approx P_{\mathrm{thr}}$ after
applying the intensity stabilisation feedback loop and terminating it with
spread in condensate node densities of $\approx 1\%\;\mathrm{(RSD)}$.
Comparison of the density spread and the spatial coherence in both polariton
lattices without ($\epsilon=0$) and with ($\epsilon>0$) is shown in Fig. 1.

##  APPENDIX C FARFIELD INTERFERENCE OF POLARITON LATTICES

The time-averaged far-field interference pattern $\tilde{I}(\mathbf{k})$ of a
set of $N$ partially coherent, narrow-bandwidth and point-like light sources
positioned in one plane can be written in the basis of spatial frequencies
$\mathbf{k}$ as (see Appendix J)

$\tilde{I}(\mathbf{k})\propto\sum_{i,j=1}^{N}\sqrt{\tilde{I_{i}}\tilde{I_{j}}}\tilde{\mu}_{ij}\mathrm{e}^{i\mathbf{k}\cdot\mathbf{d}_{ij}}.$
(11)

For each pair $\\{i,j\\}$ of point sources $\mathbf{d}_{ij}$ denotes their in-
plane spatial separation vector and $\tilde{\mu}_{ij}$ their respective
integrated complex coherence factor. For an incoherent system
($\tilde{\mu}_{ij}=0$ for $i\neq j$) the resulting homogeneous far-field
radiation pattern is the incoherent superposition of light sources with
individual intensities $\tilde{I}_{i}$. However, for non-vanishing off-
diagonal coherence elements ($|\tilde{\mu}_{ij}|>0$) the radiation pattern
$\tilde{I}(\mathbf{k})$ described in Eq. (11) becomes inhomogeneous and is
formed by the sum over discrete Fourier-components
$\exp{(i\mathbf{k}\cdot\mathbf{d}_{ij})}$ weighted by the integrated complex
coherence factor $\tilde{\mu}_{ij}$.

Hence, in analogy to beam interference measurements of coupled laser arrays
Mahler  _et al._ (2020), we can investigate the far-field emission of the
coupled polariton condensate network to gain information about the system’s
spatial coherence properties. To remove residual PL coming from polaritons
outside the condensation (gain-)centres we spatially filter the emission of
each condensation center and measure the interference of the masked condensate
emission in the far field. As schematically illustrated in Fig. 8(a) we
project the near-field PL onto a programmable aperture which allows to
selectively mask the real space field $\Psi(\mathbf{r},t)\rightarrow
A(\mathbf{r})\Psi(\mathbf{r},t)$ with amplitude function
$A(\mathbf{r})\in(0,1)$. Subsequently we image the diffraction pattern (far
field) onto a charge-coupled device (CCD) sensor. The programmable aperture
consists of a digitally controllable reflective liquid-crystal SLM as
described in reference Lousberg  _et al._ (2006). With the use of an
additional wave plate and linear polariser the apparatus transmits only
circular polarised light which is chosen to be the same as the polarisation of
the pump laser beam. By projecting the near-field PL of the polariton
condensate lattice (Fig. 8(b)) onto an aperture that transmits only the
central $\mathrm{FWHM}\approx 2\;\mathrm{\upmu m}$ of each condensation centre
(Fig. 8(c)) we record the corresponding far-field emission (Fig. 8(d)) that
consists of the interference of all condensate nodes. According to the
Fraunhofer diffraction formula the measured far-field intensity distribution
for a fully coherent wave is proportional to the squared modulus of the
Fourier transform of the optical field at the aperture location, which is
obtained from Eq. (11) for $|\tilde{\mu}_{ij}|=1$ for all pairs $\\{i,j\\}$.
Indeed, the appearance of the reciprocal triangular lattice, shown in Fig.
8(d), indicates long-range coherence across the triangular lattice of
polariton condensates.

Figure 8: Farfield interference of coherently coupled condensates. (a)
Schematic of the experimental detection setup with microcavity (MC),
microscope objective lens (MO), lenses (L1,L2), programmable aperture (PA) and
charge-coupled device (CCD). (b) Real space PL, (c) real space masking
aperture and (d) corresponding far-field interference pattern of 61
condensates arranged in a triangular geometry.

##  APPENDIX D MEASUREMENT OF COHERENCE IN POLARITON LATTICES

Figure 9: Measurement of spatial coherence in lattices of coupled polariton
condensates. By projecting the (a) real space PL of a lattice of coupled
condensates onto a programmable aperture that spatially filters the emission
to transmit only (b) condensate node $1$, (c) condensate node $2$ or (d) both
condensate nodes $1+2$ we record the respective far-field emission (e-g). The
extracted intensity profiles $\tilde{I}_{1}(k_{\parallel})$,
$\tilde{I}_{2}(k_{\parallel})$ and $\tilde{I}_{1,2}(k_{\parallel})$
perpendicular to the orientation of the interference fringes in (g) are shown
in (h) together with a fit according to Eq. (13). Pump power-dependence of the
modulus of the integrated coherence factor $|\tilde{\mu}_{12}|$ is displayed
in (i). Extracted values of the integrated complex coherence factor
$\tilde{\mu}_{1i}$ between the central condensate node $1$ and each other
condensate node $i$ are illustrated in (j) and (k) for the cases of no
condensate node density stabilisation and with node density stabilisation. All
values are extracted at a total pump power $P=1.2P_{\mathrm{thr}}$. False
colour scale shown for (a) applies to (e-g) in linear scale and to (a) in
logarithmic scale saturated below $0.01$ of the maximum count rate.

Far field measurements shown in Fig. 8 give information about the global
coherence properties of the coupled condensate network and are analogous to
time-of-flight experiments used in cold atom systems Hadzibabic  _et al._
(2006). However, precise information about the local complex coherence factor
$\mu_{ij}$ can not be obtained. The programmable aperture in our experiment
(see Fig. 8) overcomes this problem and allows for selective interference
between any pair $\left\\{i,j\right\\}$ of condensates, such that the complex
coherence factor $\mu_{ij}$ can be spatially reconstructed Lousberg  _et al._
(2006); Lundeberg  _et al._ (2007). In particular, the far-field interference
of the emission of two spatially filtered condensate nodes
$\psi_{i}(\mathbf{r},t)$ and $\psi_{j}(\mathbf{r},t)$ centred around their
respective positions $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$ yields the
measurable intensity pattern $\tilde{I}_{i+j}(\mathbf{k})$ in
$\mathbf{k}$-space as, Thompson and Wolf (1957)

$\displaystyle\tilde{I}_{i+j}(\mathbf{k})$
$\displaystyle=\tilde{I}_{i}(\mathbf{k})+\tilde{I}_{j}(\mathbf{k})$ (12)
$\displaystyle\quad+2\sqrt{\tilde{I}_{i}(\mathbf{k})\tilde{I}_{j}(\mathbf{k})}\left|\tilde{\mu}_{ij}\right|\cos\left(\mathbf{k}\cdot\mathbf{d}_{ij}+\tilde{\theta}_{ij}\right),$

where we explicitly take into account non-homogeneous intensity distributions
$\tilde{I}_{i,j}(\mathbf{k})$ due to the finite aperture sizes. In Fig. 9 we
project the triangular lattice of 61 condensates (a) onto three different
apertures filtering the emission of the central-most condensate (b), of one of
its nearest-neighbours (c) and of both condensates simultaneously (d). The
corresponding measured far-field diffraction patterns
$\tilde{I}_{1}(\mathbf{k})$, $\tilde{I}_{2}(\mathbf{k})$ and
$\tilde{I}_{1+2}(\mathbf{k})$ are shown in Figs. 9(e-g) respectively. We
extract 1D intensity profiles $\tilde{I}(k_{\parallel})$ along the direction
vector $\mathbf{k}_{\parallel}$, which is defined as co-parallel to
$\mathbf{d}_{12}$ such that
$\mathbf{k}_{\parallel}\cdot\mathbf{d}_{12}=k_{\parallel}|\mathbf{d}_{12}|$.
Equation (12) can be rewritten as,

$\displaystyle\tilde{I}_{1+2}(k_{\parallel})$
$\displaystyle=\tilde{I}_{1}(k_{\parallel})+\tilde{I}_{2}(k_{\parallel})$ (13)
$\displaystyle\quad+2\sqrt{\tilde{I}_{1}(k_{\parallel})\tilde{I}_{2}(k_{\parallel})}\left|\tilde{\mu}_{12}\right|\cos\left(k_{\parallel}|\mathbf{d}_{12}|+\tilde{\theta}_{12}\right),$

and is fitted to the experimentally extracted intensity profiles
$\tilde{I}_{1}(k_{\parallel})$, $\tilde{I}_{2}(k_{\parallel})$ and
$\tilde{I}_{1,2}(k_{\parallel})$ to yield $|\tilde{\mu}_{12}|$ and
$\tilde{\theta}_{12}$ as shown in Fig. 9(h). By extracting
$|\tilde{\mu}_{12}|$ for varying pump power we identify the threshold-like
behaviour of coherence shown in Fig. 9(i) with vanishing coherence between the
two nodes below threshold and a sharp increase at threshold. The observed
decrease of coherence for larger pump powers $P>1.2P_{\mathrm{thr}}$ is
attributed to effects involving dephasing due to increased particle
interactions, reservoir induced noise, enhanced proliferation of topological
defects, as well as to the emergence of multi-mode emission Ohadi  _et al._
(2016) reducing the fringe visibility in time-integrated measurements.

Repeating the interference measurement for the central condensate node $1$ and
each other condensate node $j$ we fully characterise the spatial coherence
properties of the central condensate node within the triangular lattice (see
Figs. 9(j,k) or Figs. 1(e,f)).

##  APPENDIX E FINITE VS INFINITE LATTICE

In experiment finite-size effects cannot be avoided whereas in theory they can
by applying periodic boundary conditions. In Fig. 10 we investigate the
difference of finite and periodic systems through simulation of Eq. (9). In
agreement with experiment, and previous theoretical observations Kalinin and
Berloff (2018), the pseudo-spin of $\tilde{\mu}_{ij}$ tilts as one approaches
the edges of the finite (damped boundary condition) condensate lattice (see
Fig. 10(c)). This effect is expected due to the non-zero flux of particles
through the boundaries of the lattice resulting in a phase gradient differing
from standing wave solutions of $\tilde{\theta}_{ij}=0,\pi$. When the system
has periodic boundary conditions (see Fig. 10(d)) the pseudo-spin tilt
vanishes and the lattice is characterised by a homogeneous standing wave
solution (zero net-flux along $x$ and $y$ directions). The simulation in both
cases uses a lattice constant of $a=12$ $\upmu$m and is performed on a
$16a\times 16a$ grid but we only plot the innermost $11a\times 11a$ of the
grid. Therefore, removing finite size effects results in a homogeneous power
distribution among all condensate nodes and thus the phases-differences stay
homogeneously $0$ or $\pi$.

We point out that by simulating a periodic system the coherence has dropped by
a small amount. This can be understood from the fact that each pulse
proliferates vortex solitons Neshev  _et al._ (2004) which, in a finite
system, can decay out through the boundary of the system (i.e., damped
boundary conditions destroy such defects). In the absence of such a decay
channel, these defect states can survive much longer in a periodic system
leading to a drop in coherence.

Figure 10: Numerical simulation of a square lattice of ballistically expanding
polariton condensates with damped (a,c) and periodic boundary conditions
(b,d). (a,b) Time-integrated and normalised real space densities for both
types of lattices with equal pump power for each condensate node. (c,d)
Extracted absolute values and phases of the complex degree of coherence
$\mu_{1j}$ between the most central condensate node (index $1$, located at
$x,y=0$) and every other condensate node (index $j$). The lattice constant is
$a=12$ $\upmu$m and size of simulation grid is $16a\times 16a$ with displayed
results are zoomed in on an area of $11a\times 11a$.

##  APPENDIX F EXCITATION PUMP POWER DEPENDENCE

In this section we describe the pump power dependence of the stabilised
triangular lattice of 61 ballistically coupled polariton condensates with
lattice constant $a=14.9\;\mathrm{\upmu m}$ shown in Figs. 1 and 2. For the
case of pulsed excitation the pump power $P$ is a time-averaged value and a
change of $P$ is equivalent to changing the peak amplitude of the sub-
picosecond laser pulse. The recorded emission patterns for pump power
$P=1.2P_{\mathrm{thr}}$ in near-field, far-field and energy resolved far-
farfield along the symmetry axis $k_{x}=0$ are illustrated in Figs. 11(a,c,e).
The pump power geometry was adjusted using the described iterative feedback
algorithm to stabilise the emission power of all condensates at threshold such
that the spread of measured condensate emission powers shows a minimum of
$1\%$ (RSD) for $P=P_{\mathrm{thr}}$. In Fig. 11(b) we show the spread of
condensate emission powers for varying total excitation pump power while
keeping the relative pump power between different nodes constant. An increase
of the spread above condensation threshold to about $20\%$ at
$P=1.2P_{\mathrm{thr}}$ originates from the different number of coupled
neighbouring sites between condensates in the centre and at the edge of the
lattice structure leading to different additional gain from coherent coupling
between condensate nodes.

In Fig. 11(d) we illustrate the pump power dependence of the total integrated
PL demonstrating a threshold level with non-linear increase of PL intensity
which we define as condensation threshold $P_{\mathrm{thr}}$. An increase in
pump power from condensation threshold $P_{\mathrm{thr}}$ to the operational
point $P=1.2P_{\mathrm{thr}}$ leads to a 9-fold increase of total PL power.

The spectrally resolved PL in reciprocal space along the symmetry axis
$k_{x}=0$ in Fig. 11(e) demonstrates PL emission at one dominant mode with
$\mathrm{FWHM}\approx 250\;\mathrm{\upmu eV}$. In agreement with the
threshold-like behaviour of the integrated PL we can see a sharp narrowing of
polariton linewidth at condensation threshold which is shown in Fig. 11(f) .

Figure 11: Pump-power dependence of a triangular lattice of polariton
condensates. (a) Recorded near-field photoluminescence of 61 polariton
condensates pumped at $1.2$ times condensation threshold. (b) Pump-power
dependence of the relative standard deviation (RSD) of the distribution of
integrated condensate emission. The system was stabilised to obtain minimum
RSD ($1\%$) at condensation threshold $P_{\mathrm{thr}}$. (c) Recorded far-
field photoluminescence for $P=1.2P_{\mathrm{thr}}$. (d) Pump-power dependence
of the integrated emission in reciprocal space. (e) Spectrally-resolved
momentum space emission along the symmetry axis $k_{x}=0$ for
$P=1.2P_{\mathrm{thr}}$. Lower polariton branch below threshold shown as red-
dashed dispersion curve. (f) Pump-power dependence of the spectral full width
at half maximum.

##  APPENDIX G TIME-RESOLVED SYNCHRONISATION OF TWO CONDENSATES

The time-resolved formation of coherence in polariton condensates under non-
resonant pulsed excitation has been investigated for single- Nardin  _et al._
(2009); Belykh  _et al._ (2013) and two-condensate systems Christmann  _et
al._ (2014); Ohadi  _et al._ (2016). In the latter case, however, a full
description of the synchronisation process between two condensates in terms of
their complex coherence factor has not been reported. Here, we explicitly
measure the (time-resolved) complex coherence factor $\mu_{12}(t)$ of two
ballistically coupled condensation centres. Pulsed excitation of two
condensates with separation distance $d_{12}=8\;\mathrm{\upmu m}$ leads to the
time-averaged near-field photoluminescence displayed in Fig. 12(a). The pump
power $P\approx 1.2P_{\mathrm{thr}}$ is the same as for the distance-
dependence shown in Fig. 3 in the main manuscript. One bright interference
peak located in-between the two condensation centres at $x,y=0$ indicates
synchronisation with vanishing phase-difference, i.e. $\tilde{\theta}_{12}=0$.
To reveal the coherence properties between the two condensates we record their
far-field interference by spatially filtering the emission of both
condensation centres (red-dashed circles) as shown in Fig. 12(b). The modulus
of the integrated coherence factor from this time-averaged interference
pattern is extracted as $|\tilde{\mu}_{12}|=0.94$. By projecting the far-field
pattern onto the entrance slit of a streak camera (time resolution $\Delta
t\approx 2\;\mathrm{ps}$) we resolve the interference and emission of each
individual condensate node in time as shown in Figs. 12(c-e). We extract the
time-dynamics of the mutual complex degree of coherence $\mu_{12}(t)$ from the
measured interference and reference signals and illustrate its phase
$\theta_{12}(t)$ and modulus $|\mu_{12}(t)|$ in Figs. 12(f) and (g),
respectively. The time-resolved occupation $\Gamma_{11}(t)+\Gamma_{22}(t)$
with $\Gamma_{ii}(t)=\langle|\psi_{i}(t)|^{2}\rangle$ of both condensates is
shown in Fig. 12(h) and reveals pulsed polariton emission with temporal width
$\approx 18\;\mathrm{ps}$ (FWHM) and peak signal at $50\;\mathrm{ps}$ after
excitation with the pump laser. We find complete synchronisation with
$\theta_{12}=0$ and a maximum degree of coherence $|\mu_{12}|=1$ appearing
earlier than the peak PL signal at $t=41\;\mathrm{ps}$ still during the
condensate growth time. The subsequent temporal decay of mutual coherence is
slower than the decay of total condensate population with fitted $1/e$ decay
times of $80\;\mathrm{ps}$ and $13\;\mathrm{ps}$, respectively.

Figure 12: Synchronisation of two ballistically coupled polariton condensates.
(a) Near-field photoluminescence of two condensates with separation distance
$d_{12}=8\;\mathrm{\upmu m}$. (b) Recorded far-field photoluminescence when
filtering the emission of the condensate centres marked with red dashed
circles in (a). (c-e) Time-resolved and normalised far-field photoluminescence
for both condensates interfering ($1+2$) and individually ($1,2$) recorded by
projecting the photoluminescence onto the entrance slit of a streak camera
(red-dashed rectangle in (b)). Time-dependencies of extracted phase difference
$\theta_{12}$, the modulus of the complex degree of coherence $|\mu_{12}|$ and
the integrated emission signal $\Gamma_{11}+\Gamma_{22}$ of both condensate
nodes are shown in (f-h). The origin for the time axis in (f-h) is defined by
the laser arrival time. The full width at half maximum of the emission signal
($\approx 18\;\mathrm{ps}$) is highlighted in light green.

##  APPENDIX H COMPARISON TO CONTINUOUS WAVE EXCITATION

We compare the results of mutual coherence between two ballistically expanding
polariton condensates pumped under sub-picosecond excitation shown in Fig. 3
to the case of pumping using continuous wave (cw) monomode laser excitation.
To prevent heating of the microcavity sample under cw excitation, the
excitation laser is modulated by an acousto-optic modulator to generate square
wave packets of $5\;\mathrm{\upmu s}$ duration at a frequency of
$10\;\mathrm{kHz}$. In Figs. 13(a) and (b) we illustrate the recorded near-
field and far-field photoluminescence of two condensate with separation
distance $d=12\;\mathrm{\upmu m}$. The system is pumped at
$P=1.6P_{\mathrm{thr}}$ with a threshold pump power $P_{\mathrm{thr}}\approx
10\;\mathrm{mW}$ per condensate node as shown in Fig. 13(c). We note that the
near-field and far-field excitation pattern in (a) and (b) are different to
the emission patterns shown in Fig. 3(a) and (b) under pulsed excitation
despite a small difference in pump spot separation distance of less than
$1\;\mathrm{\upmu m}$. In both cases the emission patterns indicate
synchronisation in a state with $\pi$ phase-difference between the two
condensate nodes. However, the outflow wavevector (or alternatively in-plane
momentum) of polaritons under pulsed excitation is larger than under cw
excitation leading to a smaller interference fringe periodicity in both real-
and momentum-space.

The extracted decay of the complex coherence factor $\mu_{12}$ with increasing
condensate separation distance under cw excitation is illustrated in Fig.
13(d). In the same plot we show the dependency of the integrated complex
coherence factor $\tilde{\mu}_{12}$ on the separation distance $d_{12}$ of two
condensates under pulsed excitation condition (as presented in Fig. 3). While
both systems demonstrate similar coherence properties for condensate
separation distances $d_{12}\leq 20\;\mathrm{\upmu m}$, the stationary
system’s coherence $|\mu_{12}|$ is enhanced for distances larger than
$20\;\mathrm{\upmu m}$ as compared to the nonstationary system under sub-
picosecond pulsed excitaiton. A Gaussian fit (see Eq. 2) of the decay of
coherence $|\mu_{12}(d_{12})|$ with increasing condensate spacing under cw
excitation yields an effective coherence length $L_{C}=40\;\mathrm{\upmu m}$.
We argue that the build-up of coherence $|\mu_{12}(t)|$ between the two
ballistically coupled condensates under pulsed excitation (as described in
Appendix G) is fast enough for small distances $d_{12}\leq 20\;\mathrm{\upmu
m}$ as to synchronise and reach the same coherence $|\mu_{12}|$ as in steady-
state operation. For larger distances, however, the increased time-of-flight
of particles travelling in-between the two condensate nodes Töpfer  _et al._
(2020) becomes noticeable and reduces the coherence factor
$|\tilde{\mu}_{12}|$ of the nonstationary system with finite life-time.

Figure 13: Ballistic coupling of two polariton condensates pumped under
continuous wave (cw) excitation. Recorded (a) near-field and (b) far-field
photoluminescence for two condensates separated by $d_{12}=12\;\mathrm{\upmu
m}$ . (c) Extracted mutual coherence factor $|\mu_{12}|$ between the two
condensates as a function of excitation pump power $P_{\mathrm{i}}$ per
condensate. Vertical arrow marks the fixed pump power ($P\approx
1.6P_{\mathrm{thr}}$) for data shown in (a,b) and for the condensate
separation distance dependence shown in (d). The data points of a dyad under
pulsed excitation (grey triangles) are illustrated for comparison to the
system under cw excitation (orange circles). Curves in (d) represent Gaussian
fits. Scale bars in (a) and (b) correspond to $10\;\mathrm{\upmu m}$ and
$1\;\mathrm{\upmu m}^{-1}$, respectively.

##  APPENDIX I BOGOLIUBOV ANALYSIS

Let us consider the continuous-wave excitation regime where
$n_{c}(\mathbf{r})=P(\mathbf{r})$ is taken time-independent. We will also
consider the ideal case of no stochasticity by setting $dW=0$. Assuming that
the reservoir $n_{x}$ follows the dynamics of the condensate we can perform an
adiabatic elimination of Eq. (9b) and keeping terms to the first order in $R$
in Eq. (9a) we have,

$i\frac{\partial\Psi}{\partial
t}=\left[-\frac{\hbar\nabla^{2}}{2m}+gP(\mathbf{r})+\alpha|\Psi|^{2}+\frac{iP}{2}\left(1-\frac{|\Psi|^{2}}{n_{s}}\right)-\frac{i\gamma}{2}\right]\Psi,$
(14)

where $n_{s}=\Gamma_{R}/R$. Performing the standard Bogoliubov treatment where
the condensate wavefunction $\Psi$ is expanded around a fixed point solution
of Eq. (14) $\Psi_{0}$ with energy $\hbar\nu$ we write,

$\Psi=\Psi_{0}(\mathbf{r})e^{-i\nu
t}\left[1+\sum_{\mathbf{q}}\left(u_{n,\mathbf{q}}e^{i\mathbf{q}\cdot\mathbf{r}+\lambda_{n}t}+v_{n,\mathbf{q}}e^{-i\mathbf{q}\cdot\mathbf{r}+\lambda_{n}^{*}t}\right)\right]$
(15)

We wish to scrutinise the dispersion of elementary excitations (Lyapunov
exponents) with complex energies $\lambda_{n}$. Substitution of Eq. (15) into
(14) and keeping only terms linear in $u_{n,\mathbf{q}}$ and
$v_{n,\mathbf{q}}$ we obtain a linearised set of equations of motion for the
disturbances. We assume that the potential $V$ and condensate $\Psi_{0}$ are
infinite and periodic such that $P(\mathbf{r})=P(\mathbf{r}+\mathbf{a})$ and
$|\Psi_{0}(\mathbf{r})|^{2}=|\Psi_{0}(\mathbf{r}+\mathbf{a})|^{2}$, where
$\mathbf{a}=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 disturbances 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})=u_{n,\mathbf{q}}(\mathbf{r}+\mathbf{a})$ and
$v_{n,\mathbf{q}}(\mathbf{r})=v_{n,\mathbf{q}}(\mathbf{r}+\mathbf{a})$.

$\mathbf{B}_{n,\mathbf{q}}(\mathbf{r})=\begin{pmatrix}u_{n,\mathbf{q}}(\mathbf{r})\\\
v_{n,\mathbf{q}}(\mathbf{r})\end{pmatrix}.$ (16)

By Fourier transforming periodic terms into the basis of reciprocal lattice
vectors we can easily solve the energies $\lambda_{n}$ belonging to
$u_{n,\mathbf{q}}$ and $v_{n,\mathbf{q}}$,

$\mathcal{L}(\mathbf{q},\Psi_{0})\mathbf{B}_{n,\mathbf{q}}=\lambda_{n}(\mathbf{q})\mathbf{B}_{n,\mathbf{q}},$
(17)

where $\mathcal{L}$ is our Bogoliubov (Lyapunov) matrix in the crystal
momentum representation,

$\displaystyle\mathcal{L}(\mathbf{q},\Psi_{0})$
$\displaystyle=\bigg{\\{}\left[\frac{\hbar}{2m}\left[\left(q_{x}-i\frac{\partial}{\partial
x}\right)^{2}+\left(q_{y}-i\frac{\partial}{\partial
y}\right)^{2}\right]+gP(\mathbf{r})-\nu+2\alpha|\Psi_{0}|^{2}\right]\hat{\sigma}_{3}$
$\displaystyle\frac{i}{2}(P-\gamma)\hat{\sigma}_{0}+\alpha\begin{pmatrix}0&\Psi_{0}^{2}\\\
-(\Psi_{0}^{*})^{2}&0\end{pmatrix}-\frac{iP}{2n_{s}}\left[2|\Psi|^{2}\hat{\sigma}_{0}+\begin{pmatrix}0&\Psi_{0}^{2}\\\
(\Psi_{0}^{*})^{2}&0\end{pmatrix}\right]\bigg{\\}}.$ (18)

The solution $\Psi_{0}$ satisfying Eq. (14) can be obtained numerically using
periodic boundary conditions. Plugging the obtained solution $\Psi_{0}$, which
corresponds to observed condensate patterns in experiment, into the eigenvalue
problem for $\mathbf{B}_{n,\mathbf{q}}$ we can finally diagonalise our system.

##  APPENDIX J FAR-FIELD DIFFRACTION FOR NARROW-BANDWIDTH PARTIALLY COHERENT
LIGHT

Let us assume the narrow bandwith optical field
$\Psi(\mathbf{r},t)=\psi(\mathbf{r,t})\exp{(-i2\pi ct/\bar{\lambda})}$ with
mean wavelength $\bar{\lambda}$ is being truncated by a thin aperture with
transmittance function $P(\mathbf{r})$ in the plane $\mathcal{A}$ such that
the field directly after the aperture is given by the product
$P(\mathbf{r})\Psi(\mathbf{r},t)$. The signal’s bandwidth $\Delta\nu$ is
assumed to be much smaller than the central frequency
$\bar{\nu}=c/\bar{\lambda}$ such that the complex amplitude
$\psi(\mathbf{r},t)$ is a slowly varying envelope in time. The average
intensity distribution $\tilde{I}(\mathbf{q})$ of the resulting far-field
diffraction pattern described in coordinate basis $\mathbf{q}$ \- realised at
the back focal plane of a thin Fourier-transforming lens with focal length $f$
\- in analogy to Schell’s theorem Goodman (2015) can be approximated as

$\displaystyle\tilde{I}(\mathbf{q})=$ (19)
$\displaystyle\iint\displaylimits_{\mathcal{A}}\iint\displaylimits_{\mathcal{A}}P^{*}(\mathbf{r}_{1})P(\mathbf{r}_{2})\sqrt{\tilde{I}(\mathbf{r}_{1})\tilde{I}(\mathbf{r}_{2})}\tilde{\mu}(\mathbf{r}_{1},\mathbf{r}_{2})\frac{\mathrm{e}^{i\frac{2\pi}{\bar{\lambda}f}\mathbf{q}\cdot\mathbf{d}_{12}}}{(\bar{\lambda}f)^{2}}\mathrm{d}\mathbf{r}_{1}\mathrm{d}\mathbf{r}_{2},$

with distance vector $\mathbf{d}_{12}=\mathbf{r}_{1}-\mathbf{r}_{2}$,
integrated complex coherence factor
$\tilde{\mu}(\mathbf{r}_{1},\mathbf{r}_{2})$ and average intensity
$\tilde{I}(\mathbf{r}_{1,2})$ of the optical field $\psi(\mathbf{r_{1,2},t})$
at locations $\mathbf{r}_{1,2}$ in the input plane $\mathcal{A}$. We further
approximate a distribution of point-like holes in the aperture, i.e.

$P(\mathbf{r})=A_{0}\sum_{i}\delta\left(\mathbf{r}-\mathbf{r}_{i}\right),$
(20)

where $A_{0}$ corresponds to the finite physical size of each hole, and we
transform into the basis of spatial frequencies
$\mathbf{k}(\mathbf{q})=2\pi\mathbf{q}/\bar{\lambda}f$ yielding

$\tilde{I}(\mathbf{k})=\frac{A_{0}^{2}}{(\lambda
f)^{2}}\sum_{i,j}\sqrt{\tilde{I}(\mathbf{r}_{i})\tilde{I}(\mathbf{r}_{j})}\tilde{\mu}(\mathbf{r}_{i},\mathbf{r}_{j})\mathrm{e}^{i\mathbf{k}\cdot\mathbf{d}_{ij}}.$
(21)

Under the assumption of a fully coherent field with constant coherence factor
$\tilde{\mu}(\mathbf{r}_{i},\mathbf{r}_{j})=1$ the resultant far-field
diffraction pattern (Eq. (21)) reduces to the well-known Fraunhofer
diffraction formula,

$\tilde{I}(\mathbf{k})=\left|\frac{A_{0}}{\lambda
f}\sum_{i}\sqrt{\tilde{I}(\mathbf{r}_{i})}\mathrm{e}^{i\mathbf{k}\cdot\mathbf{r}_{i}}\right|^{2}.$
(22)

It is apparent from Eq. (22) that in case of a periodic arrangement (lattice)
of aperture holes at locations $\mathbf{r}_{i}$ the resultant intensity
distribution $\tilde{I}(\mathbf{k})$ is formed by the squared modulus of the
discrete Fourier-transform of the optical field sampled at locations
$\mathbf{r}_{i}$.

## References

  * Pikovsky  _et al._ (2003) A. Pikovsky, J. Kurths, M. Rosenblum, and J. Kurths,  _Synchronization: a universal concept in nonlinear sciences_ , Vol. 12 (Cambridge university press, 2003).
  * Strogatz  _et al._ (2005) S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, and E. Ott, Nature 438, 43 (2005).
  * Gray  _et al._ (1989) C. M. Gray, P. König, A. K. Engel, and W. Singer, Nature 338, 334 (1989).
  * Davis  _et al._ (1995) K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995).
  * Bennett  _et al._ (2002) M. Bennett, M. F. Schatz, H. Rockwood, and K. Wiesenfeld, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 458, 563 (2002).
  * Soriano  _et al._ (2013) M. C. Soriano, J. García-Ojalvo, C. R. Mirasso, and I. Fischer, Rev. Mod. Phys. 85, 421 (2013).
  * Watts and Strogatz (1998) D. J. Watts and S. H. Strogatz,  393, 3 (1998).
  * Strogatz (2001) S. H. Strogatz, nature 410, 268 (2001).
  * Motter  _et al._ (2013) A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, Nature Physics 9, 191 (2013).
  * Romera  _et al._ (2018) M. Romera, P. Talatchian, S. Tsunegi, F. A. Araujo, V. Cros, P. Bortolotti, J. Trastoy, K. Yakushiji, A. Fukushima, H. Kubota, _et al._ , Nature 563, 230 (2018).
  * Kao  _et al._ (2016) T.-Y. Kao, J. L. Reno, and Q. Hu, Nature Photonics 10, 541 (2016).
  * Morsch and Oberthaler (2006) O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006).
  * Struck  _et al._ (2011) J. Struck, C. Ölschläger, R. Le Targat, P. Soltan-Panahi, A. Eckardt, M. Lewenstein, P. Windpassinger, and K. Sengstock, Science 333, 996 (2011).
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbène, Nature Physics 8, 267 (2012).
  * Dung  _et al._ (2017) D. Dung, C. Kurtscheid, T. Damm, J. Schmitt, F. Vewinger, M. Weitz, and J. Klaers, Nature Photonics 11, 565 (2017).
  * Amo and Bloch (2016) A. Amo and J. Bloch, Comptes Rendus Physique 17, 934 (2016).
  * Berloff  _et al._ (2017) N. G. Berloff, M. Silva, K. Kalinin, A. Askitopoulos, J. D. Töpfer, P. Cilibrizzi, W. Langbein, and P. G. Lagoudakis, Nature Materials 16, 1120 (2017).
  * Ballarini  _et al._ (2020) D. Ballarini, A. Gianfrate, R. Panico, A. Opala, S. Ghosh, L. Dominici, V. Ardizzone, M. De Giorgi, G. Lerario, G. Gigli, T. C. H. Liew, M. Matuszewski, and D. Sanvitto, Nano Letters (2020), 10.1021/acs.nanolett.0c00435.
  * Hadzibabic  _et al._ (2006) Z. Hadzibabic, P. Krüger, M. Cheneau, B. Battelier, and J. Dalibard, Nature 441, 1118 (2006).
  * Kasprzak  _et al._ (2006) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymanska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Nature 443, 409 (2006).
  * Deng  _et al._ (2007) H. Deng, G. S. Solomon, R. Hey, K. H. Ploog, and Y. Yamamoto, Phys. Rev. Lett. 99, 126403 (2007).
  * Krizhanovskii  _et al._ (2009) D. N. Krizhanovskii, K. G. Lagoudakis, M. Wouters, B. Pietka, R. A. Bradley, K. Guda, D. M. Whittaker, M. S. Skolnick, B. Deveaud-Plédran, M. Richard, R. André, and L. S. Dang, Phys. Rev. B 80, 045317 (2009).
  * Manni  _et al._ (2011) F. Manni, K. G. Lagoudakis, B. Pietka, L. Fontanesi, M. Wouters, V. Savona, R. André, and B. Deveaud-Plédran, Phys. Rev. Lett. 106, 176401 (2011).
  * Roumpos  _et al._ (2012) G. Roumpos, M. Lohse, W. H. Nitsche, J. Keeling, M. H. Szymańska, P. B. Littlewood, A. Löffler, S. Höfling, L. Worschech, A. Forchel, and Y. Yamamoto, Proceedings of the National Academy of Sciences 109, 6467 (2012).
  * Wertz  _et al._ (2010) E. Wertz, L. Ferrier, D. D. Solnyshkov, R. Johne, D. Sanvitto, A. Lemaître, I. Sagnes, R. Grousson, A. V. Kavokin, P. Senellart, G. Malpuech, and J. Bloch, Nature Physics 6, 860 (2010).
  * Caputo  _et al._ (2018) D. Caputo, D. Ballarini, G. Dagvadorj, C. Sánchez Muñoz, M. De Giorgi, L. Dominici, K. West, L. N. Pfeiffer, G. Gigli, F. P. Laussy, M. H. Szymańska, and D. Sanvitto, Nature Materials 17, 145 (2018).
  * Ohadi  _et al._ (2018) H. Ohadi, Y. del Valle-Inclan Redondo, A. J. Ramsay, Z. Hatzopoulos, T. C. H. Liew, P. R. Eastham, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. B 97, 195109 (2018).
  * Askitopoulos  _et al._ (2019) A. Askitopoulos, L. Pickup, S. Alyatkin, A. Zasedatelev, K. G. Lagoudakis, W. Langbein, and P. G. Lagoudakis, arXiv e-prints , arXiv:1911.08981 (2019), arXiv:1911.08981 [cond-mat.quant-gas] .
  * Baboux  _et al._ (2018) F. Baboux, D. D. Bernardis, V. Goblot, V. N. Gladilin, C. Gomez, E. Galopin, L. L. Gratiet, A. Lemaître, I. Sagnes, I. Carusotto, M. Wouters, A. Amo, and J. Bloch, Optica 5, 1163 (2018).
  * Richard  _et al._ (2005) M. Richard, J. Kasprzak, R. Romestain, R. André, and L. S. Dang, Phys. Rev. Lett. 94, 187401 (2005).
  * Töpfer  _et al._ (2020) J. Töpfer, H. Sigurdsson, L. Pickup, and P. Lagoudakis, Communications Physics 3, 1 (2020).
  * Tosi  _et al._ (2012) G. Tosi, G. Christmann, N. G. Berloff, P. Tsotsis, T. Gao, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Nature Physics 8, 190 (2012).
  * Cristofolini  _et al._ (2013) P. Cristofolini, A. Dreismann, G. Christmann, G. Franchetti, N. G. Berloff, P. Tsotsis, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. Lett. 110, 186403 (2013).
  * Ohadi  _et al._ (2016) H. Ohadi, R. L. Gregory, T. Freegarde, Y. G. Rubo, A. V. Kavokin, N. G. Berloff, and P. G. Lagoudakis, Phys. Rev. X 6, 031032 (2016).
  * Alyatkin  _et al._ (2020) S. Alyatkin, J. D. Töpfer, A. Askitopoulos, H. Sigurdsson, and P. G. Lagoudakis, Phys. Rev. Lett. 124, 207402 (2020).
  * Gerchberg and Saxton (1972) R. Gerchberg and W. Saxton, Optik 35, 227 (1972).
  * Kalinin and Berloff (2018) K. P. Kalinin and N. G. Berloff, New Journal of Physics 20, 113023 (2018).
  * Wouters  _et al._ (2008) M. Wouters, I. Carusotto, and C. Ciuti, Phys. Rev. B 77, 115340 (2008).
  * Wouters and Savona (2009) M. Wouters and V. Savona, Phys. Rev. B 79, 165302 (2009).
  * Szymańska  _et al._ (2006) M. H. Szymańska, J. Keeling, and P. B. Littlewood, Phys. Rev. Lett. 96, 230602 (2006).
  * Golshani  _et al._ (1999) A. Golshani, H. Pier, E. Kapon, and M. Moser, Journal of Applied Physics 85, 2454 (1999).
  * Pier  _et al._ (2000) H. Pier, E. Kapon, and M. Moser, Nature 407, 880 (2000).
  * Zhang  _et al._ (2015) M. Zhang, S. Shah, J. Cardenas, and M. Lipson, Phys. Rev. Lett. 115, 163902 (2015).
  * Fridman  _et al._ (2010) M. Fridman, M. Nixon, N. Davidson, and A. A. Friesem, Opt. Lett. 35, 1434 (2010).
  * Wouters and Carusotto (2007) M. Wouters and I. Carusotto, Phys. Rev. Lett. 99, 140402 (2007).
  * Chiocchetta and Carusotto (2013) A. Chiocchetta and I. Carusotto, EPL (Europhysics Letters) 102, 67007 (2013).
  * Pier and Kapon (1997) H. Pier and E. Kapon, Opt. Lett. 22, 546 (1997).
  * Cilibrizzi  _et al._ (2014) P. Cilibrizzi, A. Askitopoulos, M. Silva, F. Bastiman, E. Clarke, J. M. Zajac, W. Langbein, and P. G. Lagoudakis, Appl. Phys. Lett. 105, 191118 (2014).
  * Singer  _et al._ (2012) A. Singer, F. Sorgenfrei, A. P. Mancuso, N. Gerasimova, O. M. Yefanov, J. Gulden, T. Gorniak, T. Senkbeil, A. Sakdinawat, Y. Liu, D. Attwood, S. Dziarzhytski, D. D. Mai, R. Treusch, E. Weckert, T. Salditt, A. Rosenhahn, W. Wurth, and I. A. Vartanyants, Opt. Express 20, 17480 (2012).
  * Gao  _et al._ (2015) T. Gao, C. Antón, T. C. H. Liew, M. D. Martín, Z. Hatzopoulos, L. Viña, P. S. Eldridge, and P. G. Savvidis, Applied Physics Letters 107, 011106 (2015).
  * Nogrette  _et al._ (2014) F. Nogrette, H. Labuhn, S. Ravets, D. Barredo, L. Béguin, A. Vernier, T. Lahaye, and A. Browaeys, Phys. Rev. X 4, 021034 (2014).
  * Tamura  _et al._ (2016) H. Tamura, T. Unakami, J. He, Y. Miyamoto, and K. Nakagawa, Optics express 24, 8132 (2016).
  * Bobrovska  _et al._ (2014) N. Bobrovska, E. A. Ostrovskaya, and M. Matuszewski, Phys. Rev. B 90, 205304 (2014).
  * Ostrovskaya  _et al._ (2013) E. A. Ostrovskaya, J. Abdullaev, M. D. Fraser, A. S. Desyatnikov, and Y. S. Kivshar, Phys. Rev. Lett. 110, 170407 (2013).
  * Mahler  _et al._ (2020) S. Mahler, M. L. Goh, C. Tradonsky, A. A. Friesem, and N. Davidson, Phys. Rev. Lett. 124, 133901 (2020).
  * Lousberg  _et al._ (2006) G. P. Lousberg, L. D. Lundeberg, D. L. Boiko, and E. Kapon, Opt. Lett. 31, 990 (2006).
  * Lundeberg  _et al._ (2007) L. D. A. Lundeberg, G. P. Lousberg, D. L. Boiko, and E. Kapon, Applied Physics Letters 90, 021103 (2007).
  * Thompson and Wolf (1957) B. J. Thompson and E. Wolf, J. Opt. Soc. Am. 47, 895 (1957).
  * Neshev  _et al._ (2004) D. N. Neshev, T. J. Alexander, E. A. Ostrovskaya, Y. S. Kivshar, H. Martin, I. Makasyuk, and Z. Chen, Phys. Rev. Lett. 92, 123903 (2004).
  * Nardin  _et al._ (2009) G. Nardin, K. G. Lagoudakis, M. Wouters, M. Richard, A. Baas, R. André, L. S. Dang, B. Pietka, and B. Deveaud-Plédran, Phys. Rev. Lett. 103, 256402 (2009).
  * Belykh  _et al._ (2013) V. V. Belykh, N. N. Sibeldin, V. D. Kulakovskii, M. M. Glazov, M. A. Semina, C. Schneider, S. Höfling, M. Kamp, and A. Forchel, Phys. Rev. Lett. 110, 137402 (2013).
  * Christmann  _et al._ (2014) G. Christmann, G. Tosi, N. G. Berloff, P. Tsotsis, P. S. Eldridge, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, New Journal of Physics 16, 103039 (2014).
  * Goodman (2015) J. W. Goodman,  _Statistical optics_ (John Wiley & Sons, 2015).