Theory of hydrodynamic phenomena in optical mesh lattices

Signatures of superfluid-like behaviour have recently been observed experimentally in a nonlinear optical mesh lattice, where the arrival time of optical pulses propagating in a pair of coupled optical fiber loops is interpreted as a synthetic spatial dimension. Here, we develop a general theory of the fluid of light in such optical mesh lattices. On the one hand, this theory provides a solid framework for an analytical and numerical interpretation of the experimental observations. On the other hand it anticipates new physical effects stemming from the specific spatio-temporally periodic geometry of our set-up. Our work opens the way towards the full exploitation of optical mesh lattices system as a promising platform for studies of hydrodynamics phenomena in fluids of light in novel configurations.


I. INTRODUCTION
In recent years, there has been significant growth in the field of "fluids of light", which aims to bring ideas from quantum gases into nonlinear optics [1].Although historically based on analogies between paraxial light propagation in nonlinear systems and the Gross-Pitaevskii equation (GPE) of weakly-interacting atomic gases [2][3][4][5][6], experiments in this field expanded dramatically in the mid-2000s thanks to the observation of exciton-polariton condensates in semiconductor microcavities [7].This led to rapid progress in observing superfluidity effects and topological excitations of exciton-polariton fluids [8][9][10], and has fuelled wider interest in hydrodynamic effects in light [11].This interest has also grown even further over the last decade, thanks to developments in artificial magnetic fields for light [12] and cavity-less propagating geometries [13][14][15][16].
Building on previous experiments with a focus on quantum-optical effects [17,18], we have recently shown that fluids of light can also be studied in a so-called optical mesh lattice [19].The central idea of this set-up is to exploit the arrival time of classical optical pulses propagating in coupled optical fiber loops to encode one (or more [20,21]) discrete synthetic spatial dimensions [17,18,22].In contrast to more traditional platforms for fluids of light [1,[14][15][16], the flexibility of our set-up allows for the observation of the full light-field dynamics for long propagation times; the engineering of arbitrary dynamical potentials, such as both stationary and moving defects; and the control of the effective nonlinearity in a wide range with standard optoelectronic components.Previously, these advantages have been exploited in studies of, e.g., P T -symmetric physics [18,[21][22][23][24][25] and topological effects [26][27][28][29].In our recent experiment [19], we observed several key qualitative signatures of superfluid-like behaviour, including a non-zero speed of sound due to nonlinear effects and the breakdown of apparent superfluidity above a critical-velocity threshold.
In the present paper, we develop a general and com-plete theory for fluids of light in optical mesh lattices and we present its application to our experimental set-up.
The full development of this theory requires a detailed study of a number of features that are specific to the novel platform, namely the effect of the peculiar connectivity of the optical mesh lattice and the Floquet nature of the evolution given by the repeated propagation of the light pulses around the fiber loops.In addition to providing a solid interpretation of the experimental observations in [19], we describe new superfluid hydrodynamics effects which will directly follow from the new features of our platform.These include a more complex dependence of the speed of sound on the density, the onset of dynamical instabilities of the superfluid at high densities and a weakened superfluidity by Umklapp processes due to the spatial periodicity.The structure of the article is built with two types of readers in mind.Firstly, researchers interested in reproducing the experiment and/or understanding all details of its theoretical interpretation will find that the full article is a reference summarizing all important aspects of optical mesh lattices for hydrodynamic experiments.Secondly, readers that are interested in the general physics of fluids of light in optical mesh lattices can focus on the general theory and discussion of analogies and differences with other platforms for fluids of light, but may skip those subsections that are mostly devoted to the details of realistic experimental set-ups.Specifically, in Sec.II, we present an in-depth analysis of the evolution equations for the nonlinear optical mesh lattice.Building atop the known theory, we present a complete theory of photonic bands in the linear regime and of the Bogoliubov collective excitations in a weakly interacting fluid of light.Analogies with relativistic Dirac-type superfluids are also drawn.In Sec.III, we numerically demonstrate how to accurately measure the corresponding speed of sound from the effects of a stationary defect in a realistic optical mesh lattice experiment.While Subsec.III A is of general interest, the following Subsecs.III B-III C address all those details that are of crucial interest for experiments.This provides a comprehensive theoretical framework to the results previously presented in the main text and the supplemental material of Ref. [19].In Sec.IV, we numerically investigate the effective critical velocity for superfluidity and, in particular, the impact of Umklapp processes due to the periodicity of the optical mesh lattice.To this purpose, we explore a restricted Landau criterion and we identify the apparent threshold for the emission of collective excitations by a moving defect.Once again, the first subsections Subsec.IV A-IV B are of general interest, with specific experimental details being discussed in Subsec.IV C. Finally, we draw conclusions in Sec.V.

II. NONLINEAR OPTICAL MESH LATTICES
In this section, we shall begin by introducing light evolution in nonlinear optical mesh lattices.We then review the derivation of the Bogoliubov dispersion for weak perturbations propagating on top of a stationary and uniform optical field.As we show, this predicts a speed of sound for low-wavelength and low-energy excitations, which tends to zero in the limit of vanishing nonlinearity.As we discuss, this system shares many features with a Dirac-type relativistic superfluid, such as the existence of an instability for high defocusing nonlinearity and a maximum speed of sound.The nonlinear optical mesh lattice therefore opens the way towards the experimental investigation of nonlinear and relativistic effects.

A. Evolution in a nonlinear optical mesh lattice
The discrete time-evolution of light moving in a onedimensional optical mesh lattice is realised via a timemultiplexing scheme, based on two coupled optical fiber loops [17,18,22], as shown in Fig. 1(a).In the basic set-up, a light pulse is injected into one loop, and it then propagates around the loop until it is split by the (50/50) beamsplitter into two pulses, one circulating in each loop.When these pulses again reach the beamsplitter, they are each split into two more pulses and so on.A detector in one of the loops then observes a train of pulses over time, as shown in the example in Fig. 1(b).
The time-multiplexing scheme relies on choosing the lengths of the loops, L 1 and L 2 , such that there is a small but non-zero length difference, ∆L = L 1 −L 2 , as well as a long average length L = (L 1 +L 2 )/2, with L ≫ ∆L.This guarantees that there is a clear separation of time-scales between the average round-trip time, T = (T 1 +T 2 )/2 and the relative time-delay, ∆T = T 1 − T 2 , where T 1 (T 2 ) are the times taken for a light pulse to travel around the long (short) loop.Provided that the pulse width is narrower than the minimum pulse separation ∆T , the arrival time of the pulses from each loop (before they are combined) at the beamsplitter, can be expressed as T = m T + n∆T /2, where m and n are two integers.Physically, the integer m counts the total number of round trips, while n counts how many more round trips each pulse made in the long instead of the short loop.This identification is unambiguous provided that all the pulses from each round-trip fit in the time-window T 2 set by propagation around the shorter loop; this typically allows for evolutions over several hundred round-trips depending on the length of the fiber loops chosen.At longer times, the identification of individual pulses may not be possible as a pulse with a large positive value of n after m round-trips may overlap with a pulse with a large negative value of n after m + 1 round-trips.However, experimentally the evolution of the light-field can be confined to small |n| values through the control e.g. of additional amplitude modulators thus realizing absorbing boundary conditions, in which case it is possible to propagate for an even larger number of round-trips [19,[31][32][33].
To re-interpret the evolution in terms of lattice dynamics, we note that, as the light propagates, the integer m always increases for each successive round trip, while the integer n can either increase or decrease by one depending on whether the short or long loop is traversed.Guided by this, the integer m can be interpreted as a discrete time step, while n is a discrete position index along a "synthetic spatial dimension".This leads to the effective (1 + 1)D optical mesh lattice shown schematically in Fig. 1(c).Using this mapping, we can then replot the measured pulse intensities extracted from a time trace [Fig.1(b)], to reveal the light dynamics in the synthetic optical mesh lattice [see Fig. 1(d)].
In all experiments pulses were long enough and the loop lengths similar enough that dispersion effects do not play a significant role for the experimentally realised propagation length: the pulses emerging after each round-trip have an almost identical shape independently of their previous path along short and long loops, which guarantees their perfect overlap at the beam-splitter.Under these assumptions, each pulse can be characterized by a single complex amplitude only, as further discussed e.g. in Refs.[23,31].The pulse distribution at step m+1 is then described by the evolution equations [23,30,34]: where the two components u m n and v m n of the wavefunction denote the pulse amplitude entering the beamsplitter from the short and long loop, respectively.The relative phase shift of π/2 for light that couples from one loop to the other loop is a consequence of energy conservation at the symmetric 50/50 beamsplitter.Note that an unusual feature of the optical mesh lattice is its diamond connectivity over time, as physically, light either travels through the short loop (n → n − 1) or the long loop (n → n + 1).
To reach the nonlinear optical regime, dispersion compensating optical fibers can be used for the fiber loops, di erent nomenclature distinguishes between the propagation of a classical wave in the Light Walk and the quantum mechanical nature of e.g.single photons or atoms of the Quantum Walk [2].Throughout the dissertation optical amplifiers and brig sources are used, and thus particular quantum optical e ects as e.g.non-classical number distributions or correlation e ects are not considered.For this reason, t Light Walk is used to stress out our basically classical approach of the description o the electromagnetic wave.When starting with a single lattice site excitation, it h formally shown, that the resulting propagation is identical for both classical and q mechanical light interference [13].

Implementation of the Pyramid of Beam Splitters Time Multiplexing Setup
Instead of using beam splitter cubes separated by free space, as it was suggested in realized e.g. in [12], for this project fibre optical components are chosen [W1, 2], wher propagate from one fibre coupler to the next, guided through an optical fibre.A fused fibre couplers are used in optical telecommunication systems, and are thus ava high quantity and at low costs, the demand of components for a pyramid of beam sp di erent nomenclature distinguishes between the propagation of a classical wave in the Light Walk and the quantum mechanical nature of e.g.single photons or atoms of the Quantum Walk [2].Throughout the dissertation optical amplifiers and brigh sources are used, and thus particular quantum optical e ects as e.g.non-classical p number distributions or correlation e ects are not considered.For this reason, th Light Walk is used to stress out our basically classical approach of the description of the electromagnetic wave.When starting with a single lattice site excitation, it ha formally shown, that the resulting propagation is identical for both classical and qu mechanical light interference [13].

Implementation of the Pyramid of Beam Splitters b Time Multiplexing Setup
Instead of using beam splitter cubes separated by free space, as it was suggested in realized e.g. in [12], for this project fibre optical components are chosen [W1, 2], where propagate from one fibre coupler to the next, guided through an optical fibre.Alt fused fibre couplers are used in optical telecommunication systems, and are thus avail high quantity and at low costs, the demand of components for a pyramid of beam spli di erent nomenclature distinguishes between the propagation of a classical wave in case of the Light Walk and the quantum mechanical nature of e.g.single photons or atoms in case of the Quantum Walk [2].Throughout the dissertation optical amplifiers and bright laser sources are used, and thus particular quantum optical e ects as e.g.non-classical photon number distributions or correlation e ects are not considered.For this reason, the term Light Walk is used to stress out our basically classical approach of the description of the of the electromagnetic wave.When starting with a single lattice site excitation, it has been formally shown, that the resulting propagation is identical for both classical and quantum mechanical light interference [13].

Implementation of the Pyramid of Beam Splitters by a Time Multiplexing Setup
Instead of using beam splitter cubes separated by free space, as it was suggested in [6] and realized e.g. in [12], for this project fibre optical components are chosen [W1, 2], where pulses propagate from one fibre coupler to the next, guided through an optical fibre.Although fused fibre couplers are used in optical telecommunication systems, and are thus available in high quantity and at low costs, the demand of components for a pyramid of beam splitters is To illustrate the basic operation of this set-up, we plot here the time-dependent light intensity signal (upper blue line) measured in the short loop in a simple example of "Light Walk" experiment, where a single optical pulse is initially injected into one loop of a linear system [30,31].Over time, a train of pulses is observed and the average pulse height of each pulse is extracted (lower orange line).Amplifiers are present in both loops in order to compensate for losses during the experiment.As can be seen, a distinct group of pulses appears after each time period T , corresponding to the average round-trip time.(In this example, T ≈ 5200ns as the average loop length was roughly 1km.)When we zoom in, we see that pulses within a group are further separated by multiples of ∆T (here, ∆T ≈ 35.4ns) corresponding to the time delay for a pulse to travel around the long rather than the short loop.Thanks to a clear separation of the timescales T and ∆T , we can label each pulse by (m, n), where these are integers counting, respectively, the total number of round trips and the excess number of round trips in the long compared to the short loop.Note that, as can be seen in the inset, for the present "Light Walk" experiment there is destructive interference between the different optical paths to (m, n) = (3, 0), so that the central peak is absent.(c) A schematic showing the reinterpretation of the two integers, n and m, as, respectively, the discrete position in a 1D lattice, n, and the discrete time step, m.Under this mapping, completing a round-trip in the short or long loop in (a) corresponds respectively to travelling from northeast or northwest to the southwest or southeast in the diamond lattice sketched in (c).(d) The average pulse heights extracted from the example in panel (b) can be replotted in terms of m and n to reveal the time-evolution of the light field in the optical mesh lattice.For this example, the light spreads out with the characteristic distribution of a "Light Walk", with interference effects clearly being visible.Note that due to the intrinsic diamond connectivity, only either even or odd lattice sites are physically accessible at each time step.Panels (b) and (d) are reproduced from Ref. [19].
as these have a significantly higher nonlinearity and suppress the growth of modulational instabilities inside the pulses when the light propagates over long distances with high intensity [35].In our formalism, the coefficient Γ accounts for the effective nonlinear action accumulated within a single round trip, and is taken to be the same for both loops.Note that a positive Γ (as describes usual dispersion-compensating fibers) actually corresponds to having a negative interaction energy in the language of quantum fluids.The phase-shifts φ m n and ϕ m n are controlled by phasemodulators inserted in each loop.As these modulators can respond faster than the minimum pulse separation, ∆T , the imposed phase shifts can be designed with an arbitrary dependence on both pulse indices, m and n.For particles moving in a lattice, this is analogous to adding an effective potential, which varies both with position and time step, and which can moreover be different for the two components u m n , v m n of the wave-function.This provides a versatile way to control the properties of optical mesh lattices; for example, previous experiments have utilised specially-designed phase together with amplitude modulations as a tool to create and investigate PT -symmetric optical mesh lattices [18,24,25,36], to = 0 and linear Γ = 0 optical mesh lattice, as calculated using pictures (a)-(c) respectively.Note that (f) shows the quasi-energy in the moving frame θ ′ (for the laboratory frame, see Fig. 3).While these methods all give physically consistent results, each has its own subtleties (see main text).In brief, the varying scales of the axes indicate the different sizes of the first Floquet-Bloch Brillouin zone for each method.At the Floquet-Bloch Brillouin zone boundary, corresponding states are indicated by yellow circles and blue stars, showing significant differences between each method.Firstly, for Method 1 in (d), the two bands cross at the zone boundaries, indicating a missed symmetry.Secondly, for Method 2 in (e), the bands are non-degenerate and do not touch, but the total number of states has doubled as auxiliary degrees of freedom have been included beyond the physical ones.Thirdly, for Method 3 in (f), the upper-band (lower-band) state at Q = −π/2 is continuously connected to the lower-band (upper-band) state at Q = π/2, indicating an unusual boundary condition (see Fig. 3).
investigate optical diametric drive acceleration [23], to engineer and map out the Berry curvature of the optical band-structure [37] and to confine (i.e.trap) the optical field around a desired position, n [19].In this work, we use the phase-shifts as a way to imprint defects on the light field in Sec.III and Sec.IV.

B. Photonic band structure
In order to physically understand the behaviour of the mesh lattice, the first step of our study is to look for the band structure for vanishing nonlinearity Γ = 0 and vanishing phases φ m n = ϕ m n = 0.However, there are important subtleties in this derivation stemming from the unusual diamond connectivity of the lattice with respect to the discrete time step (Fig. 1(b)).To illustrate this, we shall now discuss, in turn, three distinct theoretical approaches, based on introducing, respectively: (1) double steps in time and space, (2) an extended lattice with auxiliary lattice sites and (3) a moving frame.Each of these approaches gives physically equivalent results, but with subtleties that need to be taken into account.

Calculation method 1: Double steps
Given the diamond connectivity shown in Fig. 1(c), the most simple approach is to anticipate that the optical mesh lattice is periodic under double steps in the discrete time step m → m + 2 and the position n → n + 2. The linear evolution equations for a double time step from (1)-(2) are: where we have set Γ = 0 and φ m n = ϕ m n = 0.One can then look for plane-wave solutions of the form [18,22]: with wavevector Q and Floquet quasi-energy ϑ.Insertion of this Floquet-Bloch ansatz into the double-step equations, leads to the characteristic equation for the band structure which can be solved to give the dispersion relation [22]: where j = ±1 is the band index and the wavevector Q and Floquet quasi-energy ϑ are defined within the Floquet-Brillouin zone −π/2 ≤ Q, ϑ < π/2, reflecting that we are considering double steps in both position and time.This dispersion is plotted in Fig. 2(d); however, as can be seen, there is a band-crossing (indicated by yellow circles) at the Floquet-Brillouin zone boundary (i.e. at |Q| = |ϑ| = π/2), which hints at some missed symmetry in this method.This remaining symmetry can be broken and the band-degeneracy is lifted if we impose explicit double-step modulations in the evolution, for example through the phases, φ m n , ϕ m n .In such cases, an even cleaner approach can be used to redefine m as the number of double round-trips, such that the light evolution is described by a set of four coupled equations (instead of two).However, as we do not consider explicit doublestep modulations in this paper, we will focus on other calculation approaches, which do not give these spurious band-crossings.

Calculation method 2: Auxiliary lattice
The second method we will discuss is to artificially extend the lattice by adding in auxiliary sites (red squares) at odd (even) positions at even (odd) time steps, as shown in Fig. 2(b).When taken together with the physical lattice sites (white circles), this means that we have a full square lattice (in terms of discrete time-position space), except with only diagonal connectivity (dotted and solid lines).Thanks to this connectivity, a physical field which is initialized at m = 0 with a non-vanishing amplitude only in the physical even sites will continue to only be non-zero on physical even (odd) sites at all even (resp.odd) time steps m.
Since the system of physical plus auxiliary sites is a square lattice, we can apply the ansatz (5) to the singlestep evolution equations ( 1)-( 2), taking Γ = 0 and φ m n = ϕ m n = 0.The corresponding characteristic equation for the band structure is which gives the same form for the dispersion as Eq. 7, except now defined within −π ≤ ϑ, Q < π as plotted in Fig. 2(e).Note that the periodicity of the dispersion in the Bloch momentum and the Floquet quasi-energy is now consistent with the unit step periodicity of the lattice in both spatial and temporal directions and there are no spurious band crossings.However, as we still find the same form for the dispersion despite having doubled the domain of 2), we have actually found twice as many states as before.This can be understood by remembering that we artificially doubled our degrees of freedom by adding in the nonphysical auxiliary lattice sites.Therefore, while the dispersion relation for this extended lattice appears to be simple, the condition that the field is non-zero only on the physical sites (while vanishing on the auxiliary sites) must be reflected in the momentum-space picture in a subtle way.In fact, it is not encoded in the band dispersion (or in the selection of a particular branch), but rather it can be expressed in terms of a restriction on the amplitude of the field in the different Bloch modes of the extended lattice.To see this, we note from Eq. 8 that the eigenmodes satisfy the symmetry condition It is then easy to verify that the momentum-space field amplitudes of a generic physical initial state must satisfy in order that the field vanishes on the auxiliary sites, which were only introduced for mathematical convenience.The relation θ −j (Q + π) = θ j (Q) + π between the Floquet quasi-energies of the two bands guarantees that the condition ( 10) remains fulfilled at arbitrary even time steps m during the time-evolution.At odd time steps, instead, this same condition guarantees that the field amplitude is non-vanishing only at odd sites (i.e. the physical lattice sites).As Eq.( 10) implies, any physical state must contain pairs of Bloch eigenstate of the extended lattice, and any generic physical state will contain several Floquet quasienergies.This pairing of states removes the unphysical auxiliary degrees of freedom, making the state-counting for physical states consistent again with that in Method 1.As an explicit example of this state-pairing, let us consider the concrete example of a field that is uniform throughout the whole physical 1D lattice (i.e. that at, an even time step, has an equal amplitude at all even positions and zero amplitude at all odd positions).In terms of the Bloch eigenstates of the extended lattice, this corresponds to the superposition of the two states Q = 0, j = 1 and Q = π, j = −1 for the upper band or of the two states Q = 0, j = −1 and Q = π, j = 1 for the lower band.
Introducing auxiliary lattice sites is therefore a useful and simple method for carrying out calculations, but care is required when interpreting the results, as we shall return to again when discussing Bogoliubov calculations in the next section.We also note that as the specific connectivity of the extended lattice keeps the physical and unphysical auxiliary sectors totally disconnected, there is no practical problem in restricting our attention to one of the two components only and keeping in mind only at the end that the physical field has to contain an implicit multiplication by mod(n + m, 2).2(c), and the colours of the curves can be seen as representing the symmetry condition noted in (9).With respect to this larger momentum range, the dispersion is periodic and has normal boundary conditions at the zone boundaries (as indicated by green diamonds and pink squares).

Calculation method 3: A moving frame
The final approach we shall discuss is to carry out the calculation in a quasi-moving frame [see Fig. 2(c)].Here we fix the position zero to multiples of the round trip time of the long loop T 1 and not to the average round-trip time T as was done before.This avoids the difficulties of the two previous methods but introduces a new subtlety.To transform into this quasi-moving frame, we redraw the lattice and change the coordinate system to (m, l) where l = n−m.In this frame, the nonlinear equations (Eqs. 1 and 2) with vanishing phases φ m n = ϕ m n = 0 become: and so we only need to consider even values of l in our lattice.Considering the linear case (Γ = 0), we can solve these equations by using the Floquet-Bloch ansatz (5), making the replacement ϑ → ϑ ′ , as we are now solving for the quasi-energy ϑ ′ in the moving frame.The characteristic equation for the band structure is then as plotted in Fig. 2(f).Note now −π/2 ≤ Q < π/2 and −π ≤ ϑ ′ < π, as our equations ( 11)-( 12) describe two steps along l and one step along m.As can be seen, the dispersion has two slanted bands as it is in a moving frame.
To transform back into the laboratory frame, we can make use of a Galilean transform ϑ = ϑ ′ + Q where the moving frame is moving at a speed of one site per time step, as can be seen from Fig. 2(c).Then from Eq. 13 we recover the dispersion (7) found in Methods 1 and 2, except now with the Bloch momentum defined over −π/2 ≤ Q < π/2 and the quasi-energy over −π ≤ ϑ < π as shown in Fig. 3(a).This confirms that Method 3 is physically consistent with other methods, but with the advantage of avoiding the issues with spurious bandcrossings and over-counting of states.
However, Method 3 does have its own important subtlety, as there are now unusual boundary conditions at the Brillouin zone boundary stemming from the Galilean transform.As can be seen comparing Fig. 2(f) and Fig. 3(a), we must identify the upper (lower) band eigenstate at Q = −π/2 with the lower (upper) band eigenstate at Q = π/2.In fact, this is the same symmetry condition as (9), and corresponds to a shift of π in the quasi-energy when we match-up opposite edges of this Brillouin zone.This becomes crucial when we use the laboratory-frame dispersion [Fig.3(a)] to construct the extended-zone dispersion in Fig. 3(b), corresponding to the first two Brillouin zones.Here, we see that by using the unusual boundary conditions we recover the same dispersion as in Method 2 [Fig.2(e)].An advantage of using Method 3 as compared to Method 2 is that the symmetry condition on the eigenstates ( 9) is immediately clear from the colouring of the bands.Note also that for this larger momentum range, −π ≤ Q < π, we recover normal boundary conditions, as shown by green diamonds and pink squares in Fig. 3(b).

C. Bogoliubov dispersion
In order to derive the Bogoliubov dispersion of collective excitations, one needs to consider the dynamics of small perturbations, δu m n and δv m n , on top of a stationary and initially unperturbed light field, i.e. such that [38,39].We shall now show how the Bogoliubov dispersion can be obtained by extending either the above Method 2 (of auxiliary lattice sites) or the above Method 3 (of a moving frame).

Bogoliubov theory from Method 2
Continuing from Sec. II B, we assume the unperturbed light field is at Q = 0, and is described by: where I 0 is the light intensity in each loop.The two signs correspond to the eigenstates of the j = ±1 upper and lower bands, and we have inserted the form of the corresponding eigenstates and quasi-energies ϑ = ±π/4 at Q=0 (see Eq. 7) [31].At Q = 0, the group velocity of the linear dispersion relation vanishes, and hence the light field is stationary, while the effective mass, m * ∝ (∂ 2 ϑ/∂Q 2 ) −1 , is maximal and of opposite (positive/negative) sign in the two (upper/lower) bands.
Note that this calculation considers a field that is nonvanishing at all sites and implicitly assumes that the physical field is obtained after multiplying by mod(n + m, 2).If one wished to restrict the field to the physical sites at all steps, one would have to develop a Bogoliubov theory around a time-dependent state involving two components at different Floquet quasi-energy, which is a much more complicated task.However, since the dynamics of the physical and auxiliary sites are decoupled also at nonlinear level, we anticipate that the Bogoliubov dispersion is not affected by the inclusion of the auxiliary sites.
In the nonlinear regime, the sign of the effective mass has important consequences on the stability of the system [40].In particular, if the nonlinearity, Γ, and the effective mass, m * , have the same sign, it corresponds to a focusing nonlinearity which destabilises the system, as for a BEC with attractive interactions [38].To avoid this instability, a defocusing nonlinearity is required; when the nonlinearity Γ > 0 (Γ < 0), this can be achieved by exciting the lower (upper) band eigenstate, for which the effective mass is negative (positive).From here on, we focus on the Γ > 0 case that is relevant to the experiments [19], and so only discuss the lower-band eigenstate in (14).
The small perturbations, δu m n and δv m n , are described by a temporally and spatially periodic ansatz [19]: where A u , A v , B u , B v are constant coefficients and where θ and k are now the "energy" and "Bloch momentum" of the Bogoliubov excitations.
Inserting this ansatz into the nonlinear evolution equations (1-2) and linearizing in the small perturbation, one gets the following matrix equation for the collective eigenmodes, T and the Bogoliubov matrix Solving this equation leads to the Bogoliubov dispersion relation: which consists of four branches, within the Floquet-Brillouin zone −π ≤ k, θ < π.In Fig. 4, these branches are plotted in black/red/green depending on the negative/positive/zero value of the Bogoliubov norm As usual, branches are organized in pairs with opposite quasi-momentum k and Floquet quasi-energy and opposite norm.As a consequence of the overall phase-rotation symmetry of the problem, a pair of positive-and negative-norm branches must go through the k = 0 and θ = 0 point according to Goldstone's theorem, and as can be seen here.However, as we are implicitly adding auxiliary lattice sites and working in the extended-lattice picture [see Fig. 2(b)], we have again included unphysical degrees of freedom in this calculation.Analogously to Sec.II B, there is a symmetry between states and quasi-energies, e.g.⃗ Φ −j,η (k+π) = ⃗ Φ j,η (k) and θ −j,η (k+π) = θ j,η (k)+π, where η is the sign of the Bogoliubov norm [such that the band labelling goes from bottom to top in Fig. 4 as (j, η) = (−1, −1), (−1, 1), (1, −1), (1, 1)].Again, we can circumvent these complications by instead carrying out the calculation in a moving frame.

Bogoliubov theory from Method 3
Now, following on from Sec. II B, we apply Bogoliubov theory to the moving frame.As above, the unperturbed light field is assumed to be at Q = 0 on the lower band and described by (14).The small perturbations, δu m l and δv m l , are also again described by ( 15) when we rewrite θ → θ ′ to represent the "energy" of the Bogoliubov excitations in the moving frame.Inserting this ansatz into the evolution equations (11)(12) and linearizing in the small perturbation, one gets the following Bogoliubov matrix which leads to which describes four slanted Bogoliubov bands, as shown in Fig. 5(a) for ΓI 0 = 0.2.Again, this dispersion can be transformed back to the laboratory frame by substituting θ = θ ′ + k, recapturing Eq. 17 and as shown in Fig. 5(b).As before, the major differences with Method 2 are that now −π/2 ≤ k < π/2 and −π ≤ θ < π, meaning that Method 3 only captures the desired physical states, but with the price of an unusual boundary condition that connects bands with opposite j (and the same Bogoliubov norm η) across the Brillouin zone boundary (as shown by symbols in Fig. 5).Enforcing this boundary condition when plotting an extended-zone scheme in Fig. 5(c) then exactly recovers the curves obtained via Method 2 in Fig. 4. The symmetry that relates eigenstates in different branches can also be clearly visualised in Fig. 5(c) through the matching colours of the curves.Within this extended-zone scheme (corresponding to the first two Brillouin zones of the physical lattice), the Bogoliubov dispersion is periodic in the usual sense, as shown by open grey symbols in Fig. 5(c).For this reason, hereafter, we shall focus on the Bogoliubov dispersion between −π ≤ k < π.
Through this extended theoretical discussion, we have therefore established the basis for a Bogoliubov description of fluids of light in optical mesh lattices; hereafter we will focus on the physical implications of these results, without further detailed discussion of the underlying subtleties.

D. Bogoliubov instabilities and speed of sound
As can be seen in Fig. 4, for finite nonlinearities, the Bogoliubov dispersion ( 17) can develop non-zero imaginary parts, indicating parameters for which the system is unstable.Firstly, this occurs when Γ < 0, corresponding to having a focusing nonlinearity for the lower-band eigenstate with negative effective mass, as already discussed above.Secondly, and more interestingly, the system is also unstable for a defocusing nonlinearity when Γ > 0.5.This is a peculiar feature of our two component model and is completely absent in the one-component GPE.It arises due to level attraction between the two opposite-norm Bogoliubov branches in the neighborhood of the two points k = ±π/2, leading to exceptional points in the Bogoliubov spectrum.We will numerically demonstrate the occurrence of this instability in Sec.IV B; this has not yet been observed experimentally as it requires a very high nonlinearity and hence a very strong light intensity.
In the intermediate region 0 ≤ ΓI 0 ≤ 0.5 of stability, the Bogoliubov dispersion is purely real; still, the dispersion is dramatically affected by the interactions.This can be seen by noting that in the linear regime (Γ = 0), straightforward trigonometric algebra shows that the positive-norm η > 0 components of the Bogoliubov dispersion ( 17) recover (modulo a shift by the Floquet quasi-energy of the unperturbed state θ = −π/4) the band dispersion (7) in the absence of nonlinearities.Hence, in the linear case, all branches have a parabolic dispersion around k = 0, π.
However, in the nonlinear regime, in the vicinity of k = 0, the upper (Goldstone) branch (i.e. with j = −1, η = 1 in the notation introduce above) acquires a linear form: indicating that the long wavelength excitations of the system are phonon-like, moving with a speed of sound: as indicated by the blue dashed lines in Fig. 4 and as shown in Fig. 6.Interestingly, this speed of sound appears to diverge for ΓI 0 → 1; however, the system is also unstable at such high nonlinearities, and so, in practice, the speed of sound is bound by its maximum value of v s → 1 that is reached for ΓI 0 → 0.5.This is still twice as large as the maximum speed which linear excitations on the empty lattice can acquire.As we now discuss, this is unlike the usual behaviour of superfluids, but instead can be viewed as a relativistic effect, stemming from the form of the evolution equations.

E. Comparison with non-relativistic and relativistic quantum fluids
We now compare the above results, firstly, with the usual non-relativistic Gross-Pitaevskii equation (GPE) describing, for example, weakly-interacing cold atomic Bose-Einstein condensates [38], and secondly, with a type of Dirac relativistic system, which is motivated from the continuum limit of the linear optical mesh lattice.
The GPE can be written as [38] where ψ is the wave function, M is the particle mass, g is the mean-field interaction parameter and we have set ℏ = 1.Expanding with respect to small perturbations to the light field as ψ = ψ 0 +δψ, the well-known Bogoliubov dispersion for a homogeneous system is found as [38,39] where I 0 = |ψ 0 | 2 is the unperturbed density, and ω and p are, respectively, the frequency and momentum of the elementary excitations.For small momenta, this dispersion can be expanded as ω(p) = v GPE p, where the speed of sound is given by as shown in Fig. 6.As in the optical mesh lattice, the system is unstable if g < 0, corresponding to attractive interactions between particles.However, unlike the optical mesh lattice, there is no second instability in this dispersion when g > 0, and v GPE simply increases as √ g without either a maximum or divergence at finite nonlinearity strength.To find a system which is more similar to the optical mesh lattice, we return first to the underlying evolution The phononlike linear dispersions (Eq.30) are marked by the blue dashed lines.The black/red/green (in greyscale resp.black/grey/light grey) color of each curve indicates the negative/positive/zero value of the Bogoliubov norm of the band, respectively.For αI0/M c 2 < 0 and αI0/M c 2 > 0.5, the curves develop a non-zero imaginary part, indicating dynamical instability.As can be seen, the low-momentum behaviour is qualitatively very similar to that of the optical mesh lattice (Fig. 4), while that at higher momenta deviates as the Dirac equation (Eq.24) describes a continuum model instead of a lattice.
equations ( 1)-( 2).In the absence of nonlinearities, Γ = 0, these equations can be written, after some algebra and relabelling, as [31]: where discretised derivatives with respect to the time step and position are recognisable on the LHS and RHS respectively.Taking the continuum limit, these equations can be written compactly as where we have taken ∆t = ∆x = 1, and where the two components of the vector (ψ, χ) T reflect the short and long loop degrees of freedom, (u, v) T .For compactness, we have introduced the Pauli matrices σ z and σ x .From the first-order spatial derivative on the RHS, this can be recognised as a type of relativistic Dirac equation for a two-component spinor wave-function.In the nonlinear regime, such a simple direct mapping is no longer possible.Nonetheless, motivated by this analogy, we consider a nonlinear version of (24) as where we have introduced a speed of light c = 1/ √ 2 and a mass quantifies the strength of the nonlinearity, whose form is chosen such that each component interacts with itself but not with the one in the other loop.Note that we did not introduce units for time and space in Eq. 24 and therefore the quantities introduced in the Dirac equation lack a unit as well.Similar equations have also been studied to describe cold atoms in a honeycomb lattice [41,42].As in Sec.II B, we begin by checking the linear (α = 0) dispersion by using a plane-wave ansatz: where P is the corresponding momentum and Ω is the frequency.This leads to Ω 2 = P 2 c 2 + M 2 c 4 , as expected for a relativistic Dirac model.As this is a continuum rather than a lattice model, none of the sub-lattice issues discussed in detail in Sec.II B are any longer relevant.The Bogoliubov dispersion of collective excitations can then be calculated, as before, by considering small perturbation on top of a stationary and initially unperturbed field, i.e. such that ψ → ψ 0 + δψ and χ → χ 0 + δχ.Similar to Sec.II C, we assume the unperturbed light field to have a finite power and to be at rest at P = 0, as described by: where the signs correspond to states with Ω(P = 0) = ∓M c 2 + αI 0 respectively.Note that in terms of the mesh lattice parameters introduced above, this corresponds to having Ω(P = 0) = ± √ 2 − √ 2ΓI 0 .As in the optical mesh lattice, our next step is to consider weak perturbations on top of a uniform field at rest in the lower band, i.e, at Ω(P = 0) = − √ 2 − √ 2ΓI 0 .The small perturbations, δψ(x, t) and δχ(x, t), are described by a temporally and spatially periodic ansatz: where A ψ , A χ , B ψ , B χ are constant coefficients and where ω and p are the frequency and momentum of the Bogoliubov excitations.Inserting this ansatz into the nonlinear Dirac equation (24) and linearizing in the small perturbations, one gets the following matrix equation for the collective eigenmodes, with ⃗ Φ D = (A ψ , A χ , B ψ , B χ ) T and the Bogoliubov matrix Solving this equation leads to the Bogoliubov dispersion relation: As in the nonlinear optical mesh lattice [c.f.Fig. 4], this has two positive and two negative branches, which are plotted in Fig. 7.These branches are plotted in black/red/green depending on the negative/positive/zero value of the Bogoliubov norm As before, branches are organized in pairs with opposite momentum p and frequency ω and opposite norm, and one pair of these branches goes through the p = 0 and ω = 0 point to satisfy Goldstone's theorem.As in the nonlinear optical mesh lattice, the Bogoliubov dispersion becomes unstable when αI 0 /M c 2 < 0 or αI 0 > M c 2 /2: the former instability stems from the focusing nonlinearity, while the latter one arises from a level attraction between pairs of opposite-norm Bogoliubov branches.The speed of sound for low momenta excitations is then as also indicated in Fig. 7 and plotted in Fig. 6.This is the same functional form as that of the speed of sound for the nonlinear optical mesh lattice (20).As in that system, the speed of sound diverges when αI 0 → M c 2 , but this divergence is cut off by the instability at αI 0 > M c 2 /2, and so instead there is a maximum speed of sound set by the speed of light: v D → c as αI 0 → M c 2 /2, as shown in Fig. 6.
In summary, we have found that the Bogoliubov dispersion for the Dirac model and the nonlinear optical mesh lattice are qualitatively similar in all major respects at low momentum.This means that the optical mesh lattice can provide an optical set-up to observe certain relativistic superfluid effects.However, differences naturally appear at high momentum as the Dirac model describes a continuum system while the optical mesh lattice is discrete; as a result, for the former, the frequency and momentum can increase without bound, while for the latter, the dispersion must be periodic (with the subtleties discussed in Sec.II C).

III. MEASURING THE SPEED OF SOUND
In the previous Section we laid down the general theory of sound propagation in fluids of light in optical mesh lattices and we have characterized the different regimes.In this Section, we discuss how the above predictions for the speed of sound can be experimentally verified by switching on and off a defect at rest, and observing how the resulting perturbations propagate in the fluid.We begin with numerical simulations for the idealised case of a temporally long and spatially wide defect in a uniform light field.We then bring this closer to experiment by investigating more realistic defect profiles, and asking what happens when such defects move in an expanding rather than a uniform light field.In so doing, we demonstrate that the underlying expansion speed has to be taken into account in order to accurately measure the local speed of sound.This provides provides a theoretical framework and further numerical support for the data analysis protocols that we applied in the main text and the supplemental material of our work [19].In the analysis of the experimental data, we were forced to focus on qualitative signatures in measurements of the speed of sound due to experimental uncertainties, such as in the overall value of the experimental nonlinearity strength.In the current work, instead, we provide a comprehensive theoretical framework that sets the stage for future more quantitative experimental measurements of the speed of sound.
To create a defect in the optical field, we use the timedependent phase-shifts in (1)-( 2) to imprint a Gaussian defect phase shift in both loops, according to where φ 0 is the defect amplitude, σ n (σ m ) is the standard deviation of the defect profile with respect to position (time) and n d (m d ) is the location of the defect peak amplitude in position (time).Throughout this section, we concentrate on the simplest case of a defect that is stationary with respect to the lattice, i.e. n d is a constant at all time steps, hereafter taken to be n d = 0 unless otherwise specified; the case of a moving defect will instead be studied in Sec.IV.

A. Idealised defects in uniform light fields
To verify the expected behaviour of the speed of sound, we begin by numerically studying the idealised case, (Top Row) Spatio-temporal colorplot of the numerical differential intensity (Eq.33) for a rapidly-varying, weak defect at rest in a uniform light field, with an increasing spatial width σn for a nonlinearity ΓI0 = 0.1 .The defect profile (Eq.31) is indicated by a solid black line, marking one standard deviation, σn and σm = 2, from the defect peak amplitude of φ0 = 0.01 at m d = 100 and n d = 0.For wide defects, the emission is dominated by low-momentum sound waves, while for narrow defects, there is significant emission of higher-momentum excitations beyond sound-waves.(Bottom Row) Corresponding plots of the momentum-space differential intensity spectrum ∆I k = ∆I mmax k (Eq.35) numerically calculated at mmax = 200 showing that as the defect width decreases, the range of significant Fourier components increases, indicating the importance of excitations at higher momenta.Note that the scale of the y-axis is chosen to highlight the spread of Fourier components but it artificially cuts off the dip in the spectrum at k = 0.
where the defect profile (31) is very wide and smoothly varying compared to the discrete position and time step, i.e. σ n , σ m ≫ 1.In this regime, all but the lowest-energy excitations are suppressed, leading to an emission dominated by long-wavelength sound waves.
To further simplify the data analysis, we also consider the simplest initial condition of a spatially-uniform stationary optical light field, with the form: corresponding to exciting the lower-band eigenstate at Q = 0. Here, I 0 is the initial intensity in each loop, and the factor of mod(n + m, 2) (with m = 0) highlights the connectivity of the optical mesh lattice, as physically only either even or odd lattice sites are occupied at any given time step (see Figs. 1 and 2).The propagation of excitations can be visualised through the normalized differential intensity: ∆I = I pert.− I unpert.2I 0 (33) where ) is the local intensity when a defect is present (absent).The numerical differential intensity, ∆I, is plotted in Fig. 8(a) as a function of m, n for three values of the nonlinearity parameter ΓI 0 = 0.1, 0.3, 0.5.A small defect amplitude, φ 0 = 0.01, and large defect widths, σ n = 10 and σ m = 50, have been chosen so as to create a wide, weak and slowly-varying intensity modulation in the uniform light field.As can be seen, turning on the defect leads to the symmetric emission of propagating intensity peaks from the defect position, while turning off the defect generates propagating intensity dips; of course, the order of the peaks and dips can be switched by changing the sign of the defect amplitude.
As the nonlinearity parameter, ΓI 0 , is increased, we observe two important and related effects in the emission pattern shown in Fig 8(a).Firstly, the speed of sound (Eq.20) increases with the nonlinearity, as evident from the increasing angle at which the sound waves propagate in these plots of time versus position.Secondly, the magnitude of the differential intensity decreases, indicating that the excitation of sound waves is suppressed at stronger nonlinearities for a stationary defect.
To calculate the speed of sound in the long-wavelength regime, we can simply track the propagation of the emitted intensity peaks (or dips) to estimate: where ∆n = n max (m 2 ) − n max (m 1 ) is the peak displacement over the time step interval ∆m = m 2 − m 1 .This is plotted in Fig 8(b), with m 1 = 800, m 2 = 1000, corresponding to late times after the defect has been turned off.This numerically-estimated speed of sound shows excellent agreement with the analytical expression from Bogoliubov theory (Eq.20).

B. Realistic defects in uniform light fields
While the above numerical results are encouraging, the parameters required to reach the long-wavelength regime are far beyond current experiments where propagation times are typically limited to a few hundred time steps and the initial optical field has a width only on the order of ten positions.As a result, experimentallyrealistic defects with σ n , σ m ≃ 1 will significantly excite both sound waves and higher-momenta Bogoliubov excitations.However, as we now show numerically, it is still possible to extract an accurate measurement of the speed of sound by tracking the propagation of excitations with the lowest group velocity.
For realistic defects, the oscillatory emission patterns can be much richer than for idealised defects, as shown, for example, in Fig. 9, where we progressively increase the spatial width of a rapidly-varying defect in a uniform light field.As can be seen, for wide defects (right), the emission is dominated by low-momentum sound waves, similar to Fig. 8, while for narrow defects (left panel), there are many more peaks and dips, moving with velocities greater than the speed of sound.To understand this, we define the differential intensity spectrum according to: where u m k and v m k are the Fourier transforms with respect to the 1D synthetic lattice spanned by n of the short-and long-loop amplitudes (respectively u m n and v m n ) at a given time step m.This is plotted in the bottom row of Fig. 9 for different defect widths.As can be seen, decreasing the defect width in space, increases the range of Fourier components excited, leading to a significant population of higher-momentum states.These perturbations with a much higher spatial frequency also propagate with a higher velocity than the long wavelength sound waves (see Fig. 4).Note that the spectrum is plotted here for The dramatic nonlinearity-dependence of the spatial shape of the excitation pattern emitted by a realistic defect is illustrated in Fig. 10.In the linear regime (left panel), the excitation pattern consists in a fan of fringes extending up to an outer "light-cone", determined by the maximum group velocity of the linear waves (dotted black line), For low nonlinearities (center-left panel), the excitation pattern is bounded by an inner "sound-cone" within which excitations are significantly suppressed (dashed black line) and whose position is set by the minimum group velocity of the low-momentum Bogoliubov modes around θ ≈ 0, corresponding to the speed of sound waves (20).The visible pattern then extends outwards from this inner limit because of the higher group velocity of some higher momentum modes.From the outside, the pattern for low nonlinearities is effectively limited by an outer "light-cone" (dotted black line), determined by the maximum group velocity: of the Bogoliubov dispersion.Note that due to the nonvanishing width and duration of the defect, excitations are not emitted from a single point in time and space.To take this into account, we plot the inner "sound-cone" in Fig. 10 starting from the bottom-right corner of the defect profile, as this approximately tracks emission from the defect at late times, and the outer "light-cone" starting from the top-right corner of the defect, as this represents emission at earlier times.The dependence of both the sound velocity and the maximum group velocity on the nonlinearity is shown in the right panel of Fig. 11.On the one hand, increasing the nonlinearity increases the speed of sound, shifting the FIG.10.Spatio-temporal colorplot of the differential intensity (Eq.33) for a rapidly-varying, weak defect at rest in a uniform light field, for different values of the nonlinearity parameter ΓI0 = 0, 0.2, 0.4, 0.5.The defect profile (Eq.31) is indicated by a solid green line, marking one standard deviation, σn = 1 and σm = 2, from the defect peak amplitude of φ0 = 0.01 at m d = 10 and n d = 0.The analytical speed of sound (Eq.20) and the maximum group velocity (Eq.37 and see Fig. 11) are marked with dashed and dotted black lines emanating from the bottom-right and top-right corners of the defect profile, respectively.At low nonlinearities, the emission is effectively bounded between the dashed and dotted lines, while at high nonlinearities, the emission is dominated by sound waves, which move with the maximum group velocity of all Bogoliubov excitations.inner "sound-cone" outwards.On the other hand, at low nonlinearities, the maximum group velocity of the Bogoliubov modes remains well-approximated by its value in the linear regime, max(v 0 g ) = 1/ √ 2, such that that the outer "light-cone" does not significantly change, and so the oscillating emission pattern is confined within a narrower region.
At high enough nonlinearities (ΓI 0 ≥ 1/3, center-right and right panels in Fig. 10), the speed of sound exceeds the maximum value of the group velocity in the linear regime as shown in Fig. 11, and the inner and outer bounds merge; in this limit, sound waves dominate, and the emission pattern resembles the one of long and wide defects shown in Fig. 8. Interestingly, the defect could also excite slowly-propagating modes either from the upper Bogoliubov band or from around k = ±π (see Fig. 4), which would appear within the inner sound-cone; however, the amplitudes of excitations within this region are hardly visible for the defect parameters considered in Fig. 10.
The more complicated oscillatory emission pattern for realistic defect profiles shown in Fig. 10 requires a more careful approach to quantitatively extract the speed of sound, as compared to the simple estimate provided by Eq. 34.Guided by the identification of the inner "sound cone" in this regime, we ask whether sound waves can be identified with the innermost large and smooth feature of the emission pattern.To ascertain this, we extract the speed of these features, and show that this indeed exactly matches the analytical speed of sound.To measure the speed of the innermost large perturbations, we fit the central region of the numerical data at each time step with the ansatz: where a 1 , b 1 , c 1 , d 1 are fitting parameters [19].This function is chosen as it can reproduce the largest features of the center of the perturbation pattern: namely two inner-most minima (maxima) followed by two maxima (minima) for a positive (negative) defect amplitude.An example fit is shown in Fig. 12(a), where the fitting region extends between the two innermost maxima.Using such fits, we numerically extract the positions of the inner-most large minima, as shown in Fig. 12(b).As can be seen, the position steadily increases over time corresponding to the excitations propagating outwards from the center.To then measure the speed, we fit the minima positions with the linear function: where a 2 and b 2 are fitting parameters.We can then associate the value of a 2 with the speed of sound as is plotted in Fig. 12(c).As can be seen, these numerical results are in excellent agreement with the analytical speed of sound (20), suggesting that we can indeed identify these large minima with sound waves.

C. Realistic defects in expanding light fields
A further complication comes from the fact that in the experiments [19] the initial light field is not spatially uniform, as we have assumed above, but is prepared with some localised profile.In the absence of any other potentials, the light field will therefore expand over time.This expansion is shown, for example, in Fig. 13(a) for ΓI 0 = 0.1, starting from an initially Gaussian profile: where σ G = 50 is the Gaussian width and the factor of mod(n, 2) again accounts for the diamond connectivity of the optical mesh lattice.As compared to the Light Walk shown in Fig. 1(b,d), it displays an overall Gaussian intensity profile and a spatially constant phase.Experimentally, it can be created through a protocol described in detail in Refs.[23,31].Despite this apparent similarity, the expanding light field does have two very important consequences for the propagation of sound waves.Firstly, the light field expansion will drag perturbations with it, artificially inflating the measured excitation speed.Secondly, in the absence of gain or loss, the local intensity of the light field will drop over time, causing the effective local speed of sound to vary.We now show that both of these effects can be accounted for by applying a numerical procedure to estimate the local expansion speed of the light field.We previously applied this procedure to analyse the experimental data in [19], as well as associated numerical simulations.As mentioned above, in that case, the experimental data could not be quantitatively compared to the analytical speed of sound due to experimental uncertainties, for example, in determining the value of the experimental nonlinearity parameter.In what follows, we provide a general theoretical framework for this physics.
As above, we first fit cuts of the differential intensity at each time step using Eq.38; this is shown, for example, for m = 100 in Fig. 13(c), where the blue points are the numerical data and the black curve is the fit.We extract the positions of the innermost intensity minima by finding the minima of the fits, with the results plotted in Fig. 13(d) as a function of time.As we want to allow for the possibility that the speed of the intensity minima varies over time, we perform a moving linear fit to the data, which gives us an estimate for the instantaneous speed of the minimum, v M shown by red data points in Fig. 13(e).Here, we have chosen the moving fit for a time step m to be taken over eleven time steps centred around m.As can be seen, the estimated speed is (i) not a constant as a function of time step m, and (ii) significantly larger than the expected analytical speed of sound for the initial nonlinearity (20), plotted in (e) with a blue dashed line.To explain these observations, we need to take into account the expansion of the light field over time.
As discussed in Ref. [19], our numerical practice protocol for estimating the expansion speed of the light field is based on a version of the 1D continuity equation: where v E is the drag speed, β is an overall linear amplification factor and , in the absence of the defect.The summation domain A is chosen as −∞ ≤ n ≤ A when A > 0 or A ≤ n ≤ ∞ when A < 0; this ensures that the integral always includes at least half of the light field to minimize the effects of noise.Physically, Eq. ( 41) balances the rate of change of the intensity within a selected region to the amount of light flowing in/out and the net gain within that region.The factor of 1/2 in the first term on the right-hand side of the equation takes into account that the lattice spacing is 2 in this system, and so that the local density of light is given by I(A)/2.
In the optical mesh lattice, time is discretized into discrete time steps and we can approximate this continuity equation as:  (39), taken over eleven time steps centred around m.There are big discrepancies between vM and the expected analytical value v 0 s of the speed of sound for the peak of the initial intensity (blue dashed line).(f) Estimate vs = vM − vE for the speed of sound (blue dots), which is now in excellent agreement with the analytical speed of sound vs(nmin) using the local unperturbed intensity at the position of the intensity minimum (solid black line).
where we have included explicitly the dependence on the time step.To determine the linear amplification factor β, we can apply this equation to the summed intensity over the entire system: Numerically, we choose to set β = 0, but experimentally, β can be non-zero due to intrinsic loss and gain in the system [19].Putting this back into (42), the numerically- extracted instantaneous drag speed, v E , at the position of the minimum, n min , can be calculated.The result for ∆m = 4 is shown as black dots in Fig. 13(e).Note that n min is obtained as the minimum of a fitting function and so can take continuous values, while the drag speed can only be calculated at the available lattice sites of the optical mesh lattice, making this an approximation of the speed.In any case, as expected, the drag speed increases over time, corresponding to the intensity minima moving outwards into faster-moving regions of the light field.
Importantly, we see from Fig. 13(e) that the drag speed v E accounts for a significant fraction of the observed speed of the intensity minimum, v M .In Fig. 13(f), this is shown explicitly by plotting the estimated speed of sound v s = v M − v E , which turns out much lower than v M and which decreases over time.This estimation for v s is also much lower than the analytically-expected result (20) for the initial nonlinearity, ΓI 0 = 0.1 [blue dashed line in (e)].This is easily explained in terms of the drop of the local intensity at the position of the defect over time, which leads to a decreasing effective speed of sound.For a quantitative analysis, we can extract the average local intensity in each loop at the position of the minimum I(n min )/2, and use this value within a sort of local density approximation to calculate a local instantaneous speed of sound, v s (n min ) via (20).The result is plotted with a solid black line in Fig. 13(f) and is indeed in excellent agreement with the numerical data for all times of experimental interest.This provides a solid numerical confirmation of the protocol used to extract the speed of sound in the recent experiment [19].
Differently from experiments, we can extend our numerical analysis out to even later times as shown in Fig. 14.There we see that marked deviations do occur between the numerical and the analytical local instantaneous speed of sound around 250 ≲ m ≲ 400.To under-stand this feature, in the right panels (A-D) we plot the corresponding spatial intensity profiles of the optical field in the absence of the defect.Here, vertical dashed lines indicate the positions ±n min of the innermost minima in the intensity profiles in the presence of the defect.From these plots, we infer that the deviations occur when the sound waves move through the shoulder of the optical field, where the intensity drops sharply as a function of position and the local density approximation that implicitly underlies v s (n min ) is no longer justified.On the other hand, agreement between these estimates is excellent in all other cases where the local intensity of the underlying optical field is sufficiently slowly varying, even after the defect has passed through the shoulder of the spatial intensity profile.
As further checks, we repeat the calculation for different values of the initial nonlinearity parameter, ΓI 0 , as shown in the top panels of Fig. 15.The numerical estimate of the local speed of sound is again very accurate except when the perturbation is crossing the sharp shoulder of the spatial intensity profile.This shoulder is reached earlier for higher nonlinearities due to the higher sound velocity, which leads to an earlier onset of the deviations.We can also vary the initial Gaussian width of the optical field, σ G , as shown in the bottom panels of Fig. 15.Decreasing the size of the initial field leads again to an earlier onset of the deviation window, as less time is needed for the defect to reach the shoulder of the spatial intensity profile.
Note that in the recent experiment [19], the Gaussian width and maximum time step were respectively σ G = 8 and m max = 110, due to the number of round-trips required to generate the initial Gaussian field before the start of the superfluidity experiment [23,31].Further numerical calculations with these experimentally-realistic parameters were analysed in the Supplemental Material FIG. 15. (Top Row) Plot of the different speeds as in Fig. 14 (a) but for two different values of the nonlinearity parameter ΓI0 = 0.05 and 0.15 (instead of 0.1).Again there is very good agreement between the numerically-estimated speed of sound (blue dots) and the analytical prediction for the local speed of sound vs(nmin) (black line).The deviation between these two again occurs over an interval where the defect hits the shoulder of the spatial intensity profile (as in Fig. 14 (b)); as the expansion rate of the optical field grows with the nonlinearity, the time at which the discrepancy occurs correspondingly decreases.(Bottom Row) The same plot for ΓI0 = 0.1 but narrower initial Gaussian profile widths σG = 30 and 10 (instead of 50).As the initial width of the optical field decreases, the defect hits the shoulder of the intensity profile at an earlier time, leading to an earlier deviation between the estimated and local speed of sound.For very narrow clouds, noise in the numerical estimates is more apparent as the defect quickly leaves the central region of high intensity and the very concept of local speed of sound becomes inaccurate.Nevertheless, a good qualitative agreement can still be observed.
of Ref. [19].In the future, it would be interesting to realize more roundtrips, which would leave more time for the creation of the initial field and thus allow for the creation of wider Gaussian light fields and hence more accurate experimental measurements of the local speed of sound.

IV. SUPERFLUIDITY PROPERTIES
The characterization of collective excitations and of their propagation speed reported in the previous Section is the natural starting point to investigate the superfluidity properties of the fluid of light in a nonlinear optical mesh lattice.In this Section, we summarize our advances in this direction, based on an extension of the Landau criterion for superfluidity [1,38,39] to our specific spatiotemporally periodic geometry.
The Landau criterion, in its simplest formulation, provides a general and intuitive way to predict the breakdown of superfluid behaviour when a fluid system is traversed by a uniformly moving impurity at sufficiently large speed.This simple criterion amounts to deriving a critical velocity: below which a weak moving perturbation is not able to generate collective excitations in the fluid, which then behaves as a frictionless superfluid.Above this critical velocity collective excitations start being generated in the fluid and dissipations sets in.
For the simple case of a three-dimensional, homogeneous BEC with short-range interactions described by the GPE (see Eq. 22), it can be shown that the critical velocity is directly given by the speed of sound: [38,39].This means that the breakdown of superfluidity is associated with the emission of long-wavelength (k → 0) sound waves, and that superfluidity vanishes altogether in the limit of a noninteracting Bose gas, where v c = v GPE = 0.In other systems, the emergence of additional features in the Bogoliubov dispersion significantly affects the superfluid response; for example, the Bogoliubov dispersions of superfluid helium or a dipolar BEC can both exhibit a so-called "roton" minimum at a finite momentum which then sets the critical velocity to a lower value than the speed of sound [38].In polariton light fields, there can also be important non-equilibrium effects, leading to the emergence of novel scattering regimes [43].
In this section, we explore the effect of a moving defect in an optical mesh lattice, firstly by considering suitable generalizations of the Landau criterion in Sec.IV A, and secondly from numerical simulations in Sec.IV B. As a main result of our study, we find that for any lattice, a strict application of the Landau criterion would predict a vanishing critical velocity, so that a moving defect would always excite the light field [19].However, in practice, this prediction is strongly modified when the spatial defect width is increased, so as to reduce scattering at higher momentum.In this regime, we observe in fact that the response of the optical field is effectively "superfluid-like": a characteristic threshold for the onset of dissipation is apparent, which depends on the nonlinearity as expected from a restricted Landau criterion, and which observably drops below the speed of sound for high enough nonlinearities.

A. Umklapp processes
The Landau criterion applies the principles of energy and momentum conservation to the Bogoliubov excita-tion spectrum of a system [38,39].For example, if we have a defect moving at a constant speed, s, (as realised e.g. by taking an m-dependent n d = s(m − m d ) in Eq. 31), then a simple way to visualise the Landau criterion is by drawing on top of the excitation spectrum a straight line passing through the origin with a slope s.In continuous space and time, the excitation processes that are allowed by energy and momentum conservation correspond to the intersections of this line with the Bogoliubov dispersion [44].In particular, the critical velocity, v c , is then read off as the speed at which the straight line first touches the Bogoliubov bands from below.If the system is a superfluid, then v c ̸ = 0, and for s < v c , there are no resonant excitations and, thus, no dissipation.
The situation is slightly more complicated in our spatially and temporally periodic optical mesh lattice configuration, as illustrated schematically in Fig. 16(a) and (b).Given the discrete spatial and temporal periodicities, one needs in fact to work in the extended-zone scheme and include replicas of the straight line according to these periodicities.As a result, a direct application of the Landau criterion leads to a striking conclusion: the critical velocity is always zero.This can be seen immediately from the above geometric argument by realising that the aforementioned straight lines always intersects with the Bogoliubov dispersion of the optical mesh lattice in some Brillouin zone, see e.g.Fig. 16(a) and (b).On this basis, one may expect that a lattice system can never respond like a true superfluid to a moving defect, as the defect always leads to dissipation.This notwithstanding, we note that such systems can still display other typical features of superfluidity, for example for what concerns the robustness of a superflow to dynamical and energetic instabilities [45][46][47][48] or the rigidity of the system under "twisted boundary conditions" [49].
Furthermore, the above argument based on the Landau criterion does not take into account the momentumselectivity of the excitation process due to the width of the defect, σ n .In practice, excitations at higher momenta k ≫ 1/σ n are suppressed in the case of a wide defect.To see this, let us consider a general lattice defect potential moving at a constant speed with the form V (n − sm, m) = β(m)α(n − sm), where α(n − sm) corresponds to the moving spatial profile of the defect and β(m) to the temporal profile.For example, a stationary defect of the form Eq. 31 with a m-independent position n d is captured by this functional form with Gaussianshaped α(n) and β(m) and s = 0. Imposing periodic boundary conditions, the discrete Fourier transform of the moving defect potential can be defined as: N qn (45) where N denotes the number of sites along n, and M that along m, with ε and q being integers that run from zero up to M −1 and N −1 respectively and which are related to a (discretized) energy and momentum as θ = (2π/M )ε and k = (2π/N )q respectively.Continuum variables are recovered in the limit that M and N tend to ∞.
To isolate the dependence on m in the second part of the equation, we can apply Poisson's summation formula where f (k) = ∞ −∞ dx f (x) e −ikx is the standard continuous Fourier transform of f (x), to the second part of the equation.This gives (47) where ᾱ(k) is the continuous Fourier transform of α(j).
Substituting this formula back into the full expression (45) where the sums have been extended to infinity to account for the N, M → ∞ limit, we obtain: (48) where β(θ) and β(θ) = ∞ −∞ dt β(t) e iθt are the discrete and continuous Fourier transforms of β(m), respectively.
The r ′ ̸ = 0 replicas separated by ∆θ = 2π stem from the effective Floquet periodicity of our configuration.The r ̸ = 0 replicas separated by ∆k = 2π correspond to the Bragg momentum of the lattice: in the time domain, the corresponding frequency side-bands can be physically understood as the motion of the defect making the effective potential periodically either overlap with some site or fall in the void between sites, with temporal frequency ∆θ = 2πs.
For slowly varying potentials σ m ≫ 1, the Fourierspace defect potential Ṽ (k, θ) is concentrated along an infinite series of parallel straight lines θ = (k + 2πr)s + 2πr ′ corresponding to the usual defect dispersion and its replicas, as displayed in Fig. 16(a) and (b) [50].For a given r, r ′ line, the defect weight will also be controlled by ᾱ(k + 2πr), which reaches a maximum when k = −2πr.This maximal weight is indicated schematically in Fig. 16(a) and (b) by increasing the size of the markers (which compose the straight lines) around k ≈ −2πr.As can be seen, the defect weight is most significant around the center of the green dashed box, and, due to the replicas, around equivalent points in the extended Floquet-Bloch zone.
As the argument, |k + 2πr|, quickly becomes large for Umklapp excitations, such components will only be visible for sufficiently narrow defects such that the Fourier transform of the defect potential has a significant weight at these large argument values.For wide-enough defects with 2πσ n ≥ 1, the Landau criterion restricted to a small momentum/frequency range at the center of the Brillouin zone −π < k, θ < π can instead be expected to provide a good description of the system behaviour.We can verify the validity of this argument numerically, by simulating the constant motion of spatially-narrow defects (taking n d = s(m−m d ) in Eq. 31) across a uniform light field.To minimise additional effects arising from turning on and off the defect, we choose a large temporal defect width σ m = 100 and a small defect amplitude φ 0 = 0.01.In order to identify the excitation of higher Bogoliubov bands, we first replot in Fig. 16(c) and (d) the dispersion from (a) and (b), but now showing only the corresponding fun-damental r = r ′ = 0 line so that we can clearly indicate, with numbers and letters respectively, the first expected intersections between the defect dispersion and the Bogoliubov bands.We then plot in Fig. 16 (e) and (f) the logarithm of the absolute value of the differential intensity spectrum (35) as a function of the momentum with respect to the effective 1D lattice; this is calculated numerically in each case for a spatial width of σ n = 0.25 (black line) and σ n = 0.75 (red line) at the late time step m = 900 with m d = 500.As can be seen, the largest peak occurs at k = 0, signifying the depletion of the light field into Bogoliubov quasi-particles and the residual excitation of low momenta sound waves.The momenta of the next largest peaks in the spectra are all in excellent agreement with the positions of the lowest intersection points highlighted in Fig. 16(c) and (d), and indicated in Fig. 16 (e) and (f) by labelled arrows.These higher Floquet-Umklapp excitations are much weaker than the low-energy excitations at k = 0, and they get suppressed when the spatial width of the defect σ n is increased.FIG. 17.The critical velocity as calculated from the Bogoliubov dispersion within a region −kc ≤ k ≤ kc.For kc = 2π or larger, vc = 0 for all stable nonlinearities, indicating that the system is not superfluid.However, for lower values of kc, a non-vanishing vc ̸ = 0 is found in the presence of a finite nonlinearity, indicating that "superfluid-like" behaviour may be observed for sufficiently wide defects.
As Umklapp excitations are therefore effectively strongly-suppressed for wide-enough defects, we may expect that restricting the Landau criterion to a limited momentum range may well-describe the system behaviour.In Fig. 17, we plot the critical velocity that would be predicted by applying the Landau criterion to the Bogoliubov dispersion within a restricted region −k c ≤ k ≤ k c .In the linear regime, the critical velocity is always vanishing independently of the cut-off k c , as superfluidity is a nonlinear phenomenon.For k c = 2π or larger (i.e. for a restricted region equal to twice the Brillouin zone shown in Fig. 4 or larger), the critical velocity also vanishes for all stable nonlinearities, indicating that the system is never a true superfluid.However, for lower values of k c , the critical velocity is non-zero in the presence of a nonlinearity, with v c → v s as k c → 0, as it happens in simple BECs.This means that the system may appear to be "superfluid-like" with a relative suppression of dissipation.
To confirm that such behaviour can be indeed observed in practice, we now numerically simulate the motion of moving defects in a nonlinear optical mesh lattice first for uniform and then for expanding light fields, comparing our results with the restricted Landau criterion.

B. Superfluid-like response of uniform light fields
We first investigate the superfluid-like response of a stationary and uniform optical field.To this purpose, we again apply a Gaussian phase defect, as discussed in Sec.III, but we now move the defect by changing the peak defect position, n d as a function of the time step m at a constant speed, s, across the uniform optical field, according to n d = s(m − m d ).We shall then measure how much the defect excites the system as a function of the defect speed.

Narrow defect
In Fig. 18, we start with the case of a spatially very narrow defect.Such a defect is able to excite higher-momentum perturbations.The numerically computed spatio-temporal differential intensity is plotted in Fig. 18(a) for three different values of the defect speed.As can be seen, the excitation pattern is very complicated.
To quantify the total intensity of the emission from the defect, it is then useful to introduce the observable: which temporally averages the spatially integrated differential intensity over the last j time steps of the propagation.
If we plot the total emission intensity (49) as a function of the speed s and the nonlinearity ΓI 0 , then a clear structure is observed as significant emission only occurs within certain ranges of defect speeds [Fig.18(b)], as would be expected for a standard superfluid.Outside this range of s, the defect cannot efficiently excite any perturbations in the first (i.e.innermost) Bogoliubov band, as all resonant excitations happen to be at too large momenta even for the narrow defect considered here: only excitations within −π ≤ k < π are significant and all emission lies within the window predicted by applying the Landau criterion over this limited momentum range.
The thickness of the orange stripe in Fig. 18(b), indicating a sizable total emission, depends on the curvature of the lowest Bogoliubov band, which is positive for small ΓI 0 , then negative at large ΓI 0 .As a guide to the eye, we have indicated the speed of sound (black dashed line); the lower and upper critical velocities for the first Bogoliubov band (red dotted line and blue dashed line respectively).As can be seen, at low nonlinearities, the speed of sound and the lower critical velocity of the first band coincide, and the emission is concentrated between this and the upper critical velocity.The relatively flat distribution of the emission intensity within the stripe is related to the velocity-independence of the friction force in onedimensional superfluids above the critical speed [51].In contrast, at higher nonlinearities, the speed of sound coincides with the upper critical velocity of the first band, and so emission from the first band is concentrated below this speed.For the sake of completeness, the lower critical velocity for the second Bogoliubov band has been shown in Fig. 18(b) as a yellow dashed line: except for a small region at high nonlinearities ΓI 0 ≲ 0.5, emission into this band is suppressed by the small value of the excitation matrix element.
Further insight on this physics can be obtained by comparing the total emission intensity shown in Fig. 18(b) to a numerical estimate for how many Bogoliubov states approximately satisfy our restricted Landau criterion as a function of the optical nonlinearity and defect speed.As discussed above, a simple way to visualise the restricted Landau criterion is to ask where the defect dispersion θ = ks intersects with the Bogoliubov bands with k < k c = π.To numerically estimate what fraction of available states may satisfy this criterion, we discretize the momentum, k, into N T equally-spaced values over the range k ∈ [0, π], and then calculate the associated set of discrete Bogoliubov energies Θ(k) for the two bands assuming 0 ≤ θ < π in Eq. 17.We then define the fraction of these states which approximately satisfy the restricted Landau criterion as: where the factor of 2 in the denominator occurs because here we consider the two lowest-energy Bogoliubov bands and where N [|Θ(k)−ks| < δθ] denotes the number of discrete Bogoliubov states for which the energy lies within a small energy window ±δθ of the defect dispersion.The small parameter δθ helps compensate for the discretization, but is also physically motivated by the finite energywidth expected for a defect switched on for a finite time [c.f.Eq. 48].
In Fig. 18(c), we plot Eq.50 for N T = 15000 and δθ = 0.001 (which is on the order of the energy-width expected for a Gaussian defect with σ m = 200), with the same marked lines as shown in Fig. 18(b).As can be seen, the largest fraction of available Bogoliubov states is around the speed of sound (black dashed line), where strong emission is also observed numerically in Fig. 18(b).Note that the fraction of states in Eq.50 also predicts other features like a significant emission around the upper critical velocity for the first Bogoliubov band (blue dashed line), which are not observed numerically in Fig. 18(b); this is because Eq.50 does not take into account other effects, such as the finite spatial width of the defect which suppresses larger-momentum excitations nor the k-dependence of the matrix elements.

Wide defect
We then move to the case of a wide and smoothlyvarying defect.The resulting spatio-temporal differential intensity is shown in Fig. 19(a) for three different defect speeds; as the defect is very wide and slowlyvarying (with the same parameters as in Fig. 8(a)), we are in the regime where the defect predominantly excites sound waves.The strength of emission is then strongest when the defect moves at the speed of sound, such that these excitations are resonant.This is clearly visible in Fig. 19(b), which shows how the total emission ( 49) is indeed peaked around the speed of sound (dashed black line).
Interestingly, by increasing the defect amplitude, we can also clearly observe the onset of the instability discussed in Sec.II C for a defocusing nonlinearity.This is shown in Figure 20 (a) and (b), where we compare the propagation of a weak and strong defect through a uniform light field at a high nonlinearity ΓI 0 = 0.45.As can be seen, both defects emit sound waves; however, in the case of the strong defect, the system becomes unstable, leading to a cascade of emitted excitations.We note that the instability does not appear to be triggered by the defect itself, but instead by the positive density bump that is present in the density modulation emitted by the defect.This suggests that the instability arises as a result of The total emission intensity (Eq.49), plotted as function of the defect velocity and nonlinearity for mmax = 1000 and j = 101.Excitations are mostly suppressed when the defect moves slower than the speed of sound (Eq.20, dashed black line).Note that little emission is also observed for defect speeds above the speed of sound, as then the defect is most resonant with higher-momenta excitations which remain relatively suppressed for wide defects [c.f.Eq. 48].
the local light intensity, I eff , increasing sufficiently to enter the unstable regime ΓI eff > 0.5, shown in Fig. 4.This is also consistent with Figure 20 (c), where we plot the emission for the strong defect in (b) as a function of the defect speed and nonlinearity.Unlike for a weak defect [Fig.19 (b)], we observe an unstable region of parameter space, which extends out from ΓI 0 = 0.5, centered along the speed of sound and which pushes down to even lower nonlinearities as the defect amplitude is further increased.

C. Superfluid-like response of expanding light fields
As introduced in Section III C, the initial light field in a realistic experiment is typically spatially-localised and then, unless confining potentials are applied, freely The instability appears at large nonlinearities, spreading from the regime, ΓI0 > 0.5, where the light field is itself unstable.Note that the scale of the colorbar is chosen to highlight both stable and unstable regions; in fact, the emission intensity for the latter case is several orders of magnitude larger than for the former, and diverges further with increasing mmax.
expands over time.In this final section, we simulate a moving defect in these more realistic experimental conditions, to show that also here signatures of a superfluidlike behaviour and a critical velocity can be observed.
In the two panels of Fig. 21, we plot the total differential intensity (49) for a realistic narrow and rapidlyvarying defect, across a Gaussian light field with initial widths σ G = 10 and σ G = 5.The defect profile is chosen to have σ n = 1, σ m = 10; the rapid switch on and off of the defect leads to significant additional excitations, which smear out the threshold.In addition, the expansion of the light field also leads to a clear downwards shift of the emission pattern towards lower speeds s as compared to Fig. 18.
This can be attributed to two effects; firstly, the dropping intensity of the light field means that the effective nonlinearity at the time when the defect is applied is always lower than the initial nonlinearity at the center of the cloud.Here, we assess the importance of this effect by plotting the analytical speed of sound based on the initial intensity as a dashed black line, and the speed of sound for the unperturbed intensity at the central position of the defect as a solid black line; as expected, the maximum of emission is in much better agreement with the latter rather than the former, showing that the effective nonlinearity is indeed lower than would be expected from initial conditions.Secondly, the defect moves across the spatially inhomogeneous light field, and so also crosses regions of lower intensities.This leads to a sizable emission also below the expected threshold speed, as is indeed observed.Once all these complications are taken into account, the main qualitative features of the total emission intensity diagram of Fig. 21 can be recognized and understood, confirming the important role of superfluid effects also in this expanding geometry.
While these conclusions are in qualitative agreement with the observations presented in the Supplemental Material of our previous experimental work Ref. 19, note that an alternative approach eventually turned out to be more effective in providing experimental evidence of a superfluid behaviour.As is presented in the main text of Ref. 19, expansion of the cloud could be blocked by applying an additional confining potential and the defect was then periodically moved through the cloud with a sinu-soidal temporal dependence.This ensured that the light intensity did not drop significantly over time and that the defect motion was restricted to the central region of the cloud.This allowed for the experimental observation of a clear threshold in the deposited energy as a function of the defect speed, which was interpreted as evidence of the breakdown of superfluid behaviour above a certain critical speed.While this observation is in qualitative agreement with a generalized Landau criterion based on the instantaneous speed and the local density, a quantitative analysis of the experiment will require more sophisticated theoretical tools, in particular to account for the non-uniform motion of the defect and the consequent effect of acceleration on the emission.A study along these lines will be the topic of future work.

V. CONCLUSIONS
In this paper, we have reported the development of a general theory of hydrodynamic effects in fluids of light in optical mesh lattices.This study had a twofold objective in mind.
On the one hand, our theory has allowed us to characterize effects, such as the behaviour of the speed of sound, the role of Umklapp processes in weakening superfluidity and the appearance of dynamical instabilities, which are qualitatively different from that found in typical superfluids of material or luminous particles in spatially continuous geometries, and are peculiar to the spatio-temporally periodic geometry of optical mesh lattices.
On the other hand, we have performed detailed numerical simulations to show that these features can be indeed accurately measured in an optical mesh lattice, while taking into account experimental complications such as the overall expansion of the light field and realistic defect profiles.In particular, this study provides a solid conceptual framework in support of the data analysis protocols applied in our recent experimental work [19], and lays the foundation for a next generation of quantitative experimental measurements.In particular, our work paves the way for further experimental exploration of more subtle superfluid hydrodynamic phenomena, going beyond our recent proof-of-principle experimental demonstration [19] and exploiting at full the novel and specific features of optical mesh lattices.For instance, our theoretical study has shown that new phenomena such as dynamical instabilities and effective critical velocities below the speed of sound may be observed in this platform.
As compared to other promising platforms for studies of superfluid hydrodynamics effects such as cold atoms [38] or exciton-polariton fluids [1], optical mesh lattices have demonstrated direct real-time and siteresolved access to the fluid observables as well as an analogous spatio-temporal resolution in the design and application of external potentials with arbitrary shapes without the need for any cumbersome experimental apparatus.These considerations, together with the re-cent demonstration of lattices with non-trivial geometrical [37] or topological properties [28,29], and of lattices with two [21] or even more spatial dimensions [12,52], suggest the promise of optical mesh lattices as an ideal candidate to investigate unprecedented regimes of superfluid hydrodynamics and of its interplay with the most subtle geometrical and topological features of complex lattices.

2 .Figure 2 . 2 .
Figure 2.2.: a,b) Comparison of a measured and simulated Light Walk, which is equiv an optical Galton board with 180 rows and about 150 beam splitters per row.Begi the third time step the propagation of the Light Walk is dominated by interference (see inset in the left panel), which are still visible in the last row after a propagatio time steps corresponding to a length of 180 km.For a better visualization of the high of the measurement, a logarithmic colour scale is chosen depicting three orders of mag Compared to a Random Walk, the Light Walk spreads much faster with a ballisti The measurement shows the symmetric intensity pattern of the short loop.The init is inserted into the long loop.

2 .Figure 2 . 2 .
Figure 2.2.: a,b) Comparison of a measured and simulated Light Walk, which is equiva an optical Galton board with 180 rows and about 150 beam splitters per row.Beginn the third time step the propagation of the Light Walk is dominated by interference pa (see inset in the left panel), which are still visible in the last row after a propagation time steps corresponding to a length of 180 km.For a better visualization of the high of the measurement, a logarithmic colour scale is chosen depicting three orders of magn Compared to a Random Walk, the Light Walk spreads much faster with a ballistic The measurement shows the symmetric intensity pattern of the short loop.The initia is inserted into the long loop.

2 .Figure 2 . 2 .
Figure 2.2.: a,b) Comparison of a measured and simulated Light Walk, which is equivalent to an optical Galton board with 180 rows and about 150 beam splitters per row.Beginning at the third time step the propagation of the Light Walk is dominated by interference patterns (see inset in the left panel), which are still visible in the last row after a propagation of 180 time steps corresponding to a length of 180 km.For a better visualization of the high quality of the measurement, a logarithmic colour scale is chosen depicting three orders of magnitudes.Compared to a Random Walk, the Light Walk spreads much faster with a ballistic speed.The measurement shows the symmetric intensity pattern of the short loop.The initial pulse is inserted into the long loop.

FIG. 1 .
FIG. 1. (a)Optical pulses propagate around a long and a short fiber loop, coupled by a 50/50 beamsplitter.(b) To illustrate the basic operation of this set-up, we plot here the time-dependent light intensity signal (upper blue line) measured in the short loop in a simple example of "Light Walk" experiment, where a single optical pulse is initially injected into one loop of a linear system[30,31].Over time, a train of pulses is observed and the average pulse height of each pulse is extracted (lower orange line).Amplifiers are present in both loops in order to compensate for losses during the experiment.As can be seen, a distinct group of pulses appears after each time period T , corresponding to the average round-trip time.(In this example, T ≈ 5200ns as the average loop length was roughly 1km.)When we zoom in, we see that pulses within a group are further separated by multiples of ∆T (here, ∆T ≈ 35.4ns) corresponding to the time delay for a pulse to travel around the long rather than the short loop.Thanks to a clear separation of the timescales T and ∆T , we can label each pulse by (m, n), where these are integers counting, respectively, the total number of round trips and the excess number of round trips in the long compared to the short loop.Note that, as can be seen in the inset, for the present "Light Walk" experiment there is destructive interference between the different optical paths to (m, n) = (3, 0), so that the central peak is absent.(c) A schematic showing the reinterpretation of the two integers, n and m, as, respectively, the discrete position in a 1D lattice, n, and the discrete time step, m.Under this mapping, completing a round-trip in the short or long loop in (a) corresponds respectively to travelling from northeast or northwest to the southwest or southeast in the diamond lattice sketched in (c).(d) The average pulse heights extracted from the example in panel (b) can be replotted in terms of m and n to reveal the time-evolution of the light field in the optical mesh lattice.For this example, the light spreads out with the characteristic distribution of a "Light Walk", with interference effects clearly being visible.Note that due to the intrinsic diamond connectivity, only either even or odd lattice sites are physically accessible at each time step.Panels (b) and (d) are reproduced from Ref.[19].

FIG. 2 .
FIG. 2. (a)The optical mesh lattice (re-drawn from Fig.1), with a red dashed box indicating double-steps in time and space, as used in Calculation Method 1.(b) The extended mesh lattice used in Calculation Method 2, with auxiliary sites (red squares) as well as physical sites (white circles).In the extended lattice, auxiliary sites and physical sites are always decoupled due to the diamond connectivity (dotted and solid lines).(c) The optical mesh lattice as viewed in a moving frame in Calculation Method 3, with l = n − m. (d)-(f) Plots of the Floquet quasi-energy dispersion in a uniform ϕ m n = φ m n = 0 and linear Γ = 0 optical mesh lattice, as calculated using pictures (a)-(c) respectively.Note that (f) shows the quasi-energy in the moving frame θ ′ (for the laboratory frame, see Fig.3).While these methods all give physically consistent results, each has its own subtleties (see main text).In brief, the varying scales of the axes indicate the different sizes of the first Floquet-Bloch Brillouin zone for each method.At the Floquet-Bloch Brillouin zone boundary, corresponding states are indicated by yellow circles and blue stars, showing significant differences between each method.Firstly, for Method 1 in (d), the two bands cross at the zone boundaries, indicating a missed symmetry.Secondly, for Method 2 in (e), the bands are non-degenerate and do not touch, but the total number of states has doubled as auxiliary degrees of freedom have been included beyond the physical ones.Thirdly, for Method 3 in (f), the upper-band (lower-band) state at Q = −π/2 is continuously connected to the lower-band (upper-band) state at Q = π/2, indicating an unusual boundary condition (see Fig.3).

FIG. 3 .
FIG. 3. (a) The Floquet quasi-energy dispersion obtained after transforming the dispersion in Fig. 2(f) for the moving lattice (Method 3) back into the laboratory frame via the Galilean transform θ = θ ′ + Q.The blue stars and yellow circles indicate the correspondence between states identified in Fig. 2(f).(b) Applying these boundary conditions, we can construct the dispersion in an extended-zone between −π ≤ Q < π with the first Brillouin zone marked by vertical grey dotted lines.This is fully equivalent to the dispersion obtained via Method 2, shown in Fig.2(c), and the colours of the curves can be seen as representing the symmetry condition noted in(9).With respect to this larger momentum range, the dispersion is periodic and has normal boundary conditions at the zone boundaries (as indicated by green diamonds and pink squares).

FIG. 4 .
FIG.4.The real (top panels) and imaginary (bottom panels) parts of the Bogoliubov dispersion for an unperturbed light field in the lower band of the nonlinear optical mesh lattice (Eq.17), for different values of the nonlinearity parameter: ΓI0 = −0.1,0, 0.2, 0.4, 0.5, 0.55 (left to right).The phonon-like linear dispersions (Eq.19) are marked by the blue dashed lines.The black/red/green (in greyscale resp.black/grey/light grey) color of each curve indicates the negative/positive/zero value of the Bogoliubov norm of the band, respectively.For ΓI0 < 0 and ΓI0 > 0.5, the bands develop a non-zero imaginary part, indicating dynamical instability.
FIG. 5. (a) The Bogoliubov quasi-energy dispersion (18) for ΓI0 = 0.2 in the moving frame.The continuous connection of the bands across the first Brillouin zone boundaries with the unusual boundary conditions of this Method (see also Fig. 3) are indicated by the symbols; bands with (j, η) = (−1, −1), (−1, 1), (1, −1), (1, 1) are colored in pink, blue, black and red respectively (as ordered from the bottom to the top at k = 0).(b) The dispersion of (a) Galilean-transformed back into the laboratory frame via θ = θ ′ + k.(c) The extended zone scheme between −π ≤ k < π constructed from (b) (with the first Brillouin zone marked by vertical grey dotted lines).This is fully equivalent to the Bogoliubov dispersion obtained via Method 2, shown in Fig. 4, with the colours of the curves now also showing the the symmetry condition on the states.With respect to this larger momentum range, the dispersion is periodic and has normal boundary conditions (as indicated by open grey symbols).

FIG. 6 .
FIG. 6.The speed of sound as calculated for (a) the optical mesh lattice, (b) the GPE equation, (c) a Dirac-type relativistic model, inspired by the optical mesh lattice.As can be seen, the speed of sound in the optical mesh lattice and the Dirac-type model share the same functional form, and are limited by the stability of the field, in contrast to the GPE result which increases without bound.

FIG. 7 .
FIG. 7. The real (top panels) and imaginary (bottom panels) parts of the Bogoliubov dispersion for the relativistic Dirac equation (Eq.24), for different values of the nonlinearity parameter: αI0/M c 2 = −0.1,0, 0.2, 0.4, 0.5, 0.55 (left to right).All other parameters M = −2 √ 2, c = 1/ √ 2, I0 = 1 and α = −Γ √2 are chosen in analogy with the optical mesh lattice.The phononlike linear dispersions (Eq.30) are marked by the blue dashed lines.The black/red/green (in greyscale resp.black/grey/light grey) color of each curve indicates the negative/positive/zero value of the Bogoliubov norm of the band, respectively.For αI0/M c 2 < 0 and αI0/M c 2 > 0.5, the curves develop a non-zero imaginary part, indicating dynamical instability.As can be seen, the low-momentum behaviour is qualitatively very similar to that of the optical mesh lattice (Fig.4), while that at higher momenta deviates as the Dirac equation (Eq.24) describes a continuum model instead of a lattice.
FIG. 8. (a) Spatio-temporal colorplot of the numerical differential intensity (Eq.33) for a smoothly-varying, weak defect at rest in a uniform light field, with different values of the nonlinearity parameter ΓI0 = 0.1, 0.3, 0.5.The defect profile (Eq.31) is indicated by a solid black line, marking one standard deviation, σn = 10 and σm = 50, from the defect peak amplitude of φ0 = 0.01 at m d = 500 and n d = 0.This smooth defect profile predominantly excites long-wavelength sound waves.(b) The numerical sound velocity (blue circles) as a function of the nonlinearity, extracted by tracking the propagation of sound waves in (a) via Eq.34 with m1 = 800 and m2 = 1000.This shows excellent agreement with the analytical speed of sound Eq. 20(dashed black line).

FIG. 9 .
FIG. 9.(Top Row) Spatio-temporal colorplot of the numerical differential intensity (Eq.33) for a rapidly-varying, weak defect at rest in a uniform light field, with an increasing spatial width σn for a nonlinearity ΓI0 = 0.1 .The defect profile (Eq.31) is indicated by a solid black line, marking one standard deviation, σn and σm = 2, from the defect peak amplitude of φ0 = 0.01 at m d = 100 and n d = 0.For wide defects, the emission is dominated by low-momentum sound waves, while for narrow defects, there is significant emission of higher-momentum excitations beyond sound-waves.(Bottom Row) Corresponding plots of the momentum-space differential intensity spectrum ∆I k = ∆I mmax

FIG. 11
FIG. 11. (a)The Bogoliubov dispersion for the optical mesh lattice (Eq.17), for the different values of the nonlinearity parameter ΓI0 = 0, 0.2, 0.4, 0.5 considered in Fig.10, as previously shown in Fig.4.The blue dashed lines indicate the phononlike linear dispersions with the speed of sound given by Eq.(19).(b) The absolute value of the group velocity, |vg| = |∂θ/∂k|, for the same values of the nonlinearity used in (a), shown with a solid blue (dark grey in greyscale) line for the two middle Bogoliubov bands (i.e.those meeting at θ = 0 in (a)) and shown with a solid purple (light grey in greyscale) line for the other two outer bands.Note that in the linear regime, for ΓI0 = 0, the blue and purple lines are identical.In each plot, the analytical speed of sound (Eq.20) and the maximum group velocity in the linear regime (max(v 0 g ) = 1/ √ 2) are marked with, respectively, dashed and dotted black lines.(c) Plot of the maximum group velocity extracted from (b) as a function of the nonlinearity parameter ΓI0 (solid blue line).The black lines indicate the speed of sound (dashed) and the maximum group velocity in the linear regime (dotted).

FIG. 12
FIG. 12. (a) Example cut of the normalised differential intensity at the final time step m = 150 for ΓI0 = 0.05.Blue points represent the numerical data and black line is the fitting function (38).Vertical dashed lines indicate the symmetric minima positions as extracted from the fitting function.The defect is chosen to have and temporal widths of σn = 2 and σm = 2 and a peak amplitude of φ0 = 0.01, and to be centered at m d = 11 and n d = 0. (b) The extracted absolute positions of the innermost large minima over time from the fits as in panel (a) for a range of different nonlinearities.(c) The numerical estimate for the speed of the sound (red crosses) as a function of the effective initial intensity as calculated from panel (b) by applying a linear fitting function (39) from m = 65 onwards.This is in excellent agreement with the analytical speed of sound (20), plotted with a solid black line.

Figure 13 (
b) then shows the numerical differential intensity resulting from the application of a realistic narrow defect, with σ n = 2, σ m = 2, φ 0 = 0.01, m d = 11 and n d = 0.As can be seen, this propagation pattern appears to be similar to that on top of a uniform light field [see e.g.Fig. 10], with excitations at large positions suppressed by the drop-off in the intensity of the cloud.

FIG. 13 .
FIG. 13.(a) Expansion of a Gaussian optical field over time, with ΓI0 = 0.1 and σG = 50.(b) Spatio-temporal colorplot of the numerical differential intensity after applying a realistic narrow defect, with σn = 2, σm = 2, φ0 = 0.01, m d = 11 and n d = 0. (c) Fit of a cut of (b) at m = 100: the blue points are the numerical data and the black curve is the fitting function (38).(d) Position of the innermost minimum of the fitting function at each time step.(e) Plot of the numericallyestimated drag speed vE (black dots), due to the expansion of the condensate, and of the estimated instantaneous speed of the intensity minimum, vM (red dots) as obtained from a moving linear fit(39), taken over eleven time steps centred around m.There are big discrepancies between vM and the expected analytical value v 0 s of the speed of sound for the peak of the initial intensity (blue dashed line).(f) Estimate vs = vM − vE for the speed of sound (blue dots), which is now in excellent agreement with the analytical speed of sound vs(nmin) using the local unperturbed intensity at the position of the intensity minimum (solid black line).

FIG. 14 .
FIG. 14. (Left panel) Results from Fig. 13(f) now plotted for longer times, with red data points indicating vM ; blue data points indicating the numerically extracted vs = vM − vE; blue dotted-dashed line indicating the speed of sound v 0s corresponding to the initial nonlinearity; the black solid line indicating the analytical prediction for the local speed of sound vs(nmin).As can be seen, there is excellent agreement between the numerical estimate and the analytical prediction (blue vs. black), except for the 250 ≲ m ≲ 400 interval.(Right panels) To understand this, we plot the intensity of the unperturbed optical field Ĩ = (|un| 2 + |vn| 2 )unpert.at different times as indicated in the left panel.The vertical dashed lines indicate the positions of the minima at each time.As can be seen, the deviations visible in the left panel appear where the wave emanating from the defect hits the shoulder of the intensity profile where the local intensity decreases sharply.

FIG. 16
FIG. 16.(a) & (b) The Bogoliubov dispersion for ΓI0 = 0.25 in an extended Floquet-Bloch zone scheme, with the range shown in Fig. 4 outlined by a green dashed box.Also plotted are the straight lines (blue points) corresponding to a defect moving at (a) s = 1.5v sound and (b) s = 0.5v sound .Due to the spatial and temporal periodicities of the optical mesh lattice, there are an infinite number of replicas of the defect straight lines.However, as discussed in the main text, for a sufficiently slowly varying defect, the corresponding Fourier-space defect potential is highly peaked around certain momenta for each line and becomes small elsewhere; this relative importance is schematically represented by an increased size of the blue points, with the largest weight being concentrated around k = 0 and equivalent points in the extended Floquet-Bloch scheme.(c) & (d) The dispersion as in (a) & (b) respectively, but only with the fundamental r = r ′ = 0 defect line (now unadjusted for the relative defect weight), in order to highlight the first intersections between these lines and the Bogoliubov bands as labelled by numbers and letters respectively.(e) & (f) The logarithm of the absolute value of the differential intensity spectra (Eq.35), as calculated at the late time step m = 900 for a defect with spatial width σn = 0.25 (black line) and σn = 0.75 (red line) moving through a uniform light field with ΓI0 = 0.25 at the speed values (e) s = 1.5v sound and (f) s = 0.5v sound .The defect parameters are: σm = 100, m d = 500 and φ0 = 0.01.Labelled arrows indicate the momenta of the lowest intersections shown in (c) & (d) back-folded into the range: −π/2 ≤ k < π/2.

s = 1 .FIG. 18
FIG. 18.(a) Spatio-temporal colorplot of the numerical differential intensity (Eq.33) for a temporally smoothly-varying, spatially narrow and very weak defect (σm = 200, σn = 1, φ0 = 0.00001, m d = 500, n d = 0), moving at speeds s = {1, 1.5, 2}v sound as indicated by the black parallelogram across a uniform light field, with nonlinearity ΓI0 = 0.1.(b) The total emission intensity (Eq.49) plotted as function of the defect velocity and nonlinearity for mmax = 1000 and j = 101.As a guide to the eye, we have marked the speed of sound (black dashed line); lower and upper critical velocities for the first Bogoliubov band (red dotted line and blue dashed line respectively); the lower critical velocity for the second Bogoliubov band (yellow dashed line) within −π ≤ k, θ < π.(c) The fraction of states (Eq.50) in the Bogoliubov bands which approximately satisfy the restricted Landau criterion as calculated with NT = 15000 and δθ = 0.001.Other marked lines are the same as in panel (b).

FIG. 19
FIG. 19.(a) Spatio-temporal colorplot of the numerical differential intensity (Eq.33) for a smoothly-varying, weak and wide defect, moving at speeds s = {0.5, 1, 1.5}v sound across a uniform light field, with nonlinearity ΓI0 = 0.25, and with all other parameters as in Fig. 8(a).The smooth defect profile (indicated by the black parallelogram) predominantly excites sound waves, with the strongest emission of radiation for s = v sound .(b)The total emission intensity (Eq.49), plotted as function of the defect velocity and nonlinearity for mmax = 1000 and j = 101.Excitations are mostly suppressed when the defect moves slower than the speed of sound (Eq.20, dashed black line).Note that little emission is also observed for defect speeds above the speed of sound, as then the defect is most resonant with higher-momenta excitations which remain relatively suppressed for wide defects [c.f.Eq. 48].

FIG. 20
FIG. 20. (a)&(b) Spatio-temporal colorplot of the numerical differential intensity (Eq.33) for, respectively, a weak defect with amplitude φ0 = 0.01, and a strong defect with amplitude φ0 = 0.2.The defect (indicated by the black parallelogram) moves at a speed s = v sound in a uniform light field, with nonlinearity ΓI0 = 0.45, and with all other parameters as in Fig. 8(a).The strong defect in (b) triggers an instability, corresponding to a cascade of emitted excitations originating from the positive density bump that is present in the excitation pattern emitted by the defect.(c) The total emission intensity (49) for the strong defect in (b), plotted as function of the defect velocity and nonlinearity for mmax = 1000 and j = 101.The instability appears at large nonlinearities, spreading from the regime, ΓI0 > 0.5, where the light field is itself unstable.Note that the scale of the colorbar is chosen to highlight both stable and unstable regions; in fact, the emission intensity for the latter case is several orders of magnitude larger than for the former, and diverges further with increasing mmax.

FIG. 21 .
FIG. 21.The total emission intensity (49) for a narrow and rapidly-varying defect, moved across a Gaussian light field with initial width σG = 10 (Top Panel) and σG = 5 (Bottom Panel).The defect parameters are σn = 1, σm = 10, m d = 20 and φ0 = 0.01.The total emission intensity (49) has been evaluated with mmax = 150 and j = 51.The analytical speed of sound from the initial intensity is shown as a dashed black line, and the speed of sound for the unperturbed intensity at the position of the defect is shown as a solid black line.As can be seen, the maximum of emission is in better agreement with the latter rather than the former, showing the impact of the light field expansion in reducing the critical speed for superfluidity.