{'fraction_non_alphanumeric': 0.036578412281739844, 'fraction_numerical': 0.01062925627481146, 'mean_word_length': 4.505982703471153, 'pattern_counts': {'":': 1, '<': 27, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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. 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 lengthL = (L 1 +L 2 )/2, withL ∆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 each pulse at the beamsplitter can be expressed as T = mT +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-windowT ; experimentally, this typically allows for evolutions over several hundred round-trips.
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 that dispersion effects could safely be neglected for the experimentally realised propagation length. Hence, each pulse is characterized by a single complex amplitude only. In that case, the pulse distribution at step m+1 is described by the evolution equations [23,30,32]: 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, 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 [33]. 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 op-  time m. Here, it is assumed that the necessary information for determining the evolution of the system is reduced to the complex field amplitudes u m n and v m n of the pulses. However, the compliance of this simplification is limited by the following practical aspects of the experiment: • The pulses possess a finite width and thus edges.
• The pulse width is related to a finite spectrum and thus subject to group velocity dispersion and higher order e ects during the propagation.
• The amplitude of each pulse is assumed to be rectangular and free of any phase structure.
The first point concerning the edge quality of the pulses can indeed influence the experiment. While the interference of the pulse plateau is in general in good accordance to the theoretical expectations based on simulating Eqs. (2.3), the edges of the pulses normally don't interfere deterministically in the experiment. However, this hurdle can in general be circumvented by averaging only over the intensity of the central part of a pulse and thus neglecting the edges. Secondly, for the linear propagation through an optical fibre, the dispersive length L D , after which a pulse of width T 0 is stretched to Ô 2T 0 , is given by [27, p. 58] di erent nomenclature distinguishes between the propagation of a classical wave in c the Light Walk and the quantum mechanical nature of e.g. single photons or atoms i 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, the 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 has formally shown, that the resulting propagation is identical for both classical and qua 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 availa high quantity and at low costs, the demand of components for a pyramid of beam split di erent nomenclature distinguishes between the propagation of a classical wave in ca the Light Walk and the quantum mechanical nature of e.g. single photons or atoms in of the Quantum Walk [2]. Throughout the dissertation optical amplifiers and bright sources are used, and thus particular quantum optical e ects as e.g. non-classical ph number distributions or correlation e ects are not considered. For this reason, the Light Walk is used to stress out our basically classical approach of the description of t the electromagnetic wave. When starting with a single lattice site excitation, it has formally shown, that the resulting propagation is identical for both classical and quan mechanical light interference [13].

Implementation of the Pyramid of Beam Splitters by Time Multiplexing Setup
Instead of using beam splitter cubes separated by free space, as it was suggested in [6] realized e.g. in [12], for this project fibre optical components are chosen [W1, 2], where p propagate from one fibre coupler to the next, guided through an optical fibre. Alth fused fibre couplers are used in optical telecommunication systems, and are thus availab high quantity and at low costs, the demand of components for a pyramid of beam splitte 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.
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 (blue line) measured in the short loop in a simple example of "Light Walk" experiment, where a single optical pulse is initially injected into a linear system [30,31]. Over time, a train of pulses is observed, with each pulse being fitted (orange line) to extract its average height. As can be seen, a distinct group of pulses appears after each time periodT , 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 timescalesT 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].
tical 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,34], to investigate optical diametric drive acceleration [23], and to engineer and map out the Berry curvature of the optical band-structure [35]. 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 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).
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(b), 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 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 Q (from −π/2 ≤ Q < π/2 in Method 1 to −π ≤ Q < π in Method 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 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).
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 timesteps, 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 con-nectivity 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).

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 timẽ 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 which leads to the dispersion 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 timestep, 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 u m n → u m n + δu m n and v m n → v m n + δv m n [36,37]. 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 [38]. 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 [36]. 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, 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 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 [36], 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 [36] i ∂ψ ∂t = − 1 2M 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 [36,37] 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 v GPE = gI 0 /M , 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 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 timestep 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 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.
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 M = −2 √ 2. The parameter α = −Γ √ 2 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 [39,40].
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: (28) 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 |A ψ | 2 + |A χ | 2 − |B ψ | 2 − |B χ | 2 . 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 further numerical verification for the data analysis protocols applied in our recent experimental work [19] and 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, where the defect profile (31)  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: where I pert. = |u m n | 2 + |v m n | 2 (or correspondingly I unpert. ) is the local intensity when a defect is present (absent). The numerical differential intensity, ∆I, is plotted in (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. 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 timesteps 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 −π/2 ≤ k < π/2 [see discussion in Sec. II C].
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 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 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 re- gion 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. 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. 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], and here we perform further numerical simulations to confirm that this method indeed accurately captures the speed of sound in an expanding light field.
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 I A = n∈A I(n) is the integrated local intensity, I(n) = |u n | 2 + |v n | 2 , 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: 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 numericallyextracted 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 understand 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 de-viations. 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.
To summarise this section, we have shown how the speed of sound can be accurately extracted using realistic defects located in expanding light fields. This provides a quantitative confirmation to the data analysis procedures recently applied to experimental data [19].

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,36,37] 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 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.
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: v c = v GPE = gI 0 /M [36,37]. 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 super-fluid 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 [36]. In polariton light fields, there can also be important non-equilibrium effects, leading to the emergence of novel scattering regimes [41].
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 excitation spectrum of a system [36,37]. 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 [42]. 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, as a moving defect always leads to dissipation. However, the above argument does not take into account the momentum-selectivity 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 Gaussian-shaped α(n) and β(m) and s = 0. Imposing periodic boundary conditions, the discrete Fourier transform of the moving defect potential can be defined as: (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ᾱ(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: 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) [43]. 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 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.
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 fundamental 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.
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 [44]. 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.  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].

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 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.
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 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. 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.
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 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 sinusoidal 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 continu-ous 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 [36] 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 recent demonstration of lattices with non-trivial geometrical [35] or topological properties [28,29], and of lattices with two [21] or even more spatial dimensions [12,45], 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.