Theory of hydrodynamic phenomena in optical mesh lattices

Signatures of superfluid-like behavior 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 spatiotemporally periodic geometry of our setup. 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.

Figure 1

(a) Optical pulses propagate around a long and a short fiber loop, coupled by a 50-50 beamsplitter. (b) To illustrate the basic operation of this setup, we plot the time-dependent light intensity signal (upper blue line) measured in the short loop in a simple example of a “Light Walk” experiment, where a single optical pulse is initially injected into one loop of a linear system [30, 31]. Over time a train of pulses is observed, and the average pulse height of each pulse is extracted (lower orange line). Amplifiers are present in both loops to compensate for losses during the experiment. As can be seen, a distinct group of pulses appears after each time period $\overline{T}$, corresponding to the average round-trip time. (In this example, $\overline{T}\approx 5200\mathrm{ns}$ as the average loop length was roughly 1 km.) When we zoom in, we see that pulses within a group are further separated by multiples of $\mathrm{\Delta }T$ (here $\mathrm{\Delta }T\approx 35.4\mathrm{ns}$) corresponding to the time delay for a pulse to travel around the long rather than the short loop. Due to a clear separation of the timescales $\overline{T}$ and $\mathrm{\Delta }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 traveling 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 (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].

Figure 2

(a) Optical mesh lattice (redrawn from Fig. 1), with a red dashed box indicating double steps in time and space, as used in Calculation Method 1. (b) The extended mesh lattice used in Calculation Method 2, with auxiliary sites (red squares) as well as physical sites (white circles). In the extended lattice, auxiliary sites and physical sites are always decoupled due to the diamond connectivity (dotted and solid lines). (c) The optical mesh lattice as viewed in a moving frame in Calculation Method 3, with $l=n-m$. (d–f) Plots of the Floquet quasi-energy dispersion in a uniform ${\varphi }_n^m={\phi }_n^m=0$ and linear $\mathrm{\Gamma }=0$ optical mesh lattice, as calculated using pictures (a)–(c), respectively. Note that (f) shows the quasi-energy in the moving frame ${\theta }^{\prime }$ (for the laboratory frame; see Fig. 3). While these methods all give physically consistent results, each has its own subtleties (see discussion). 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. First, for Method 1 in (d), the two bands cross at the zone boundaries, indicating a missed symmetry. Second, for Method 2 in (e), the bands are nondegenerate and do not touch, but the total number of states has doubled as auxiliary degrees of freedom have been included beyond the physical ones. Third, for Method 3 in (f), the upper-band (lower-band) state at $Q=-\pi /2$ is continuously connected to the lower-band (upper-band) state at $Q=\pi /2$, indicating an unusual boundary condition (see Fig. 3).

Figure 3

(a) Floquet quasi-energy dispersion obtained after transforming the dispersion in Fig. 2 for the moving lattice (Method 3) back into the laboratory frame via the Galilean transform $\theta ={\theta }^{\prime }+Q$. The blue stars and yellow circles indicate the correspondence between states identified in Fig. 2. (b) Applying these boundary conditions, we can construct the dispersion in an extended zone between $-\pi \le Q<\pi$ with the first Brillouin zone marked by vertical gray dotted lines. This is fully equivalent to the dispersion obtained via Method 2, shown in Fig. 2, and the colors 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).

Figure 4

Real (top panels) and imaginary (bottom panels) parts of the Bogoliubov dispersion for an unperturbed light field in the lower band of the nonlinear optical mesh lattice [Eq. (17)], for different values of the nonlinearity parameter: $\mathrm{\Gamma }{I}_0=-0.1,0,0.2,0.4,0.5,0.55$ (left to right). The phonon-like linear dispersions (Eq. (19) are marked by the blue dashed lines. The black, red, or green (in grayscale resp. black, gray, or light gray) of each curve indicates the negative, positive, or zero value of the Bogoliubov norm of the band, respectively. For $\mathrm{\Gamma }{I}_0<0$ and $\mathrm{\Gamma }{I}_0>0.5$, the bands develop a nonzero imaginary part, indicating dynamical instability.

Figure 5

(a) Bogoliubov quasi-energy dispersion (18) for $\mathrm{\Gamma }{I}_0=0.2$ in the moving frame. The continuous connection of the bands across the first Brillouin zone boundaries with the unusual boundary conditions of this Method (see also Fig. 3) are indicated by the symbols; bands with $(j,\eta )=(-1,-1), (-1,1), (1,-1), (1,1)$ are colored in pink, blue, black, and red, respectively (as ordered from the bottom to the top at $k=0$). (b) The dispersion of (a) Galilean transformed back into the laboratory frame via $\theta ={\theta }^{\prime }+k$. (c) The extended zone scheme between $-\pi \le k<\pi$ constructed from (b) (with the first Brillouin zone marked by vertical gray dotted lines). This is fully equivalent to the Bogoliubov dispersion obtained via Method 2, shown in Fig. 4, with the colors of the curves now also showing the symmetry condition on the states. With respect to this larger momentum range, the dispersion is periodic and has normal boundary conditions (as indicated by open gray symbols).

Figure 6

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

Figure 7

Real (top panels) and imaginary (bottom panels) parts of the Bogoliubov dispersion for the relativistic Dirac equation [Eq. (24)], for different values of the nonlinearity parameter: $\alpha I_0/M{c}^2=-0.1,0,0.2,0.4,0.5,0.55$ (left to right). All other parameters $M=-2\sqrt{2}, c=1/\sqrt{2}, I_0=1$ and $\alpha =-\mathrm{\Gamma }\sqrt{2}$ are chosen in analogy with the optical mesh lattice. The phonon-like linear dispersions [Eq. (30) are marked by the blue dashed lines. The black, red, or green (in grayscale resp. black, gray, or light gray) of each curve indicates the negative, positive, or zero value of the Bogoliubov norm of the band, respectively. For $\alpha I_0/M{c}^2<0$ and $\alpha I_0/M{c}^2>0.5$, the curves develop a nonzero imaginary part, indicating dynamical instability. As can be seen, the low-momentum behavior 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.

Figure 8

(a) Spatiotemporal colorplot of the numerical differential intensity [Eq. (33)] for a smoothly varying, weak defect at rest in a uniform light field, with different values of the nonlinearity parameter $\mathrm{\Gamma }{I}_0=0.1,0.3,0.5$. The defect profile [Eq. (31) is indicated by a solid black line, marking one standard deviation, ${\sigma }_n=10$ and ${\sigma }_m=50$, from the defect peak amplitude of ${\phi }_0=0.01$ at $m_d=500$ and $n_d=0$. This smooth defect profile predominantly excites long-wavelength sound waves. (b) The numerical sound velocity (blue circles) as a function of the nonlinearity, extracted by tracking the propagation of sound waves in (a) via Eq. (34) with $m_1=800$ and $m_2=1000$. This shows excellent agreement with the analytical speed of sound Eq. (20) (dashed black line).

Figure 9

(Top row) Spatiotemporal 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 ${\sigma }_n$ for a nonlinearity $\mathrm{\Gamma }{I}_0=0.1$. The defect profile [Eq. (31)] is indicated by a solid black line, marking one standard deviation, ${\sigma }_n$ and ${\sigma }_m=2$, from the defect peak amplitude of ${\phi }_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 $\mathrm{\Delta }{I}_k=\mathrm{\Delta }{I}_k^{m_{\text{max}}}$ [Eq. (35)] numerically calculated at $m_{\text{max}}=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$.

Figure 10

Spatiotemporal colorplot of the differential intensity [Eq. (33)] for a rapidly varying, weak defect at rest in a uniform light field, for different values of the nonlinearity parameter $\mathrm{\Gamma }{I}_0=0,0.2,0.4,0.5$. The defect profile [Eq. (31)] is indicated by a solid green line, marking one standard deviation, ${\sigma }_n=1$ and ${\sigma }_m=2$, from the defect peak amplitude of ${\phi }_0=0.01$ at $m_d=10$ and $n_d=0$. The analytical speed of sound [Eq. (20)] and the maximum group velocity [Eq. (37); see Fig. 11] are marked with dashed and dotted black lines emanating from the bottom-right and top-right corners of the defect profile, respectively. At low nonlinearities, the emission is effectively bounded between the dashed and dotted lines, while at high nonlinearities, the emission is dominated by sound waves, which move with the maximum group velocity of all Bogoliubov excitations.

Figure 11

(a) Bogoliubov dispersion for the optical mesh lattice [Eq. (17)], for the different values of the nonlinearity parameter $\mathrm{\Gamma }{I}_0=0,0.2,0.4,0.5$ considered in Fig. 10, as previously shown in Fig. 4. The blue dashed lines indicate the phonon-like linear dispersions with the speed of sound given by Eq. (19). (b) The absolute value of the group velocity, $|v_g|=|\partial \theta /\partial k|$, for the same values of the nonlinearity used in (a), shown with a solid blue (dark gray in grayscale) line for the two middle Bogoliubov bands [i.e., those meeting at $\theta =0$ in (a)] and shown with a solid purple (light gray in grayscale) line for the other two outer bands. Note that in the linear regime, for $\mathrm{\Gamma }{I}_0=0$, the blue and purple lines are identical. In each plot the analytical speed of sound [Eq. (20)] and the maximum group velocity in the linear regime ($\max \left({\mathrm{v}}_{\mathrm{g}}^0\right)=1/\sqrt{2}$) are marked with, respectively, dashed and dotted black lines. (c) Plot of the maximum group velocity extracted from (b) as a function of the nonlinearity parameter $\mathrm{\Gamma }{I}_0$ (solid blue line). The black lines indicate the speed of sound (dashed) and the maximum group velocity in the linear regime (dotted).

Figure 12

(a) Example cut of the normalized differential intensity at the final time step $m=150$ for $\mathrm{\Gamma }{I}_0=0.05$. Blue points represent the numerical data, and the black line is the fitting function (38). Vertical dashed lines indicate the symmetric minima positions as extracted from the fitting function. The defect is chosen to have spatial and temporal widths of ${\sigma }_n=2$ and ${\sigma }_m=2$ and a peak amplitude of ${\phi }_0=0.01$, and to be centered at $m_d=11$ and $n_d=0$. (b) The extracted absolute positions of the innermost large minima over time from the fits as in panel (a) for a range of different nonlinearities. (c) The numerical estimate for the speed of the sound (red crosses) as a function of the effective initial intensity as calculated from panel (b) by applying a linear fitting function (39) from $m=65$ onwards. This is in excellent agreement with the analytical speed of sound (20), plotted with a solid black line.

Figure 13

(a) Expansion of a Gaussian optical field over time, with $\mathrm{\Gamma }{I}_0=0.1$ and ${\sigma }_G=50$. (b) Spatiotemporal colorplot of the numerical differential intensity after applying a realistic narrow defect, with ${\sigma }_n=2, {\sigma }_m=2, {\phi }_0=0.01, m_d=11$, and $n_d=0$. (c) Fit of a cut of (b) at $m=100$: the blue points are the numerical data and the black curve is the fitting function (38). (d) Position of the innermost minimum of the fitting function at each time step. (e) Plot of the numerically estimated drag speed $v_E$ (black dots), due to the expansion of the condensate, and of the estimated instantaneous speed of the intensity minimum $v_M$ (red dots) as obtained from a moving linear fit (39), taken over 11 time steps centered around $m$. There are big discrepancies between $v_M$ and the expected analytical value $v_s^0$ of the speed of sound for the peak of the initial intensity (blue dashed line). (f) Estimate $v_s=v_M-v_E$ for the speed of sound (blue dots), which is now in excellent agreement with the analytical speed of sound $v_s\left(n_{\text{min}}\right)$ using the local unperturbed intensity at the position of the intensity minimum (solid black line).

Figure 14

(Left panel) Results from Fig. 13 now plotted for longer times, with red data points indicating $v_M$; blue data points indicating the numerically extracted $v_s=v_M-v_E$; blue dotted-dashed line indicating the speed of sound $v_s^0$ corresponding to the initial nonlinearity; the black solid line indicating the analytical prediction for the local speed of sound $v_s\left(n_{\text{min}}\right)$. As can be seen, there is excellent agreement between the numerical estimate and the analytical prediction (blue vs black), except for the $250\lesssim m\lesssim 400$ interval. (Right panels) To understand this, we plot the intensity of the unperturbed optical field $\tilde{I}=\left(\right|u_n{|}^2+|v_n{{|}^2)}_{\text{unpert.}}$ at different times as indicated in the left panel. The vertical dashed lines indicate the positions of the minima at each time. As can be seen, the deviations visible in the left panel appear where the wave emanating from the defect hits the shoulder of the intensity profile where the local intensity decreases sharply.

Figure 15

(Top row) Plot of the different speeds as in Fig. 14 (a) but for two different values of the nonlinearity parameter $\mathrm{\Gamma }{I}_0=0.05$ and 0.15 (instead of 0.1). Again there is very good agreement between the numerically estimated speed of sound (blue dots) and the analytical prediction for the local speed of sound $v_s\left(n_{\text{min}}\right)$ (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]; 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 $\mathrm{\Gamma }{I}_0=0.1$ but narrower initial Gaussian profile widths ${\sigma }_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.

Figure 16

(a), (b) Bogoliubov dispersion for $\mathrm{\Gamma }{I}_0=0.25$ in an extended Floquet-Bloch zone scheme, with the range shown in Fig. 4 outlined by a green dashed box. Also plotted are the straight lines (blue points) corresponding to a defect moving at (a) $s=1.5v_{\mathrm{sound}}$ and (b) $s=0.5v_{\mathrm{sound}}$. Due to the spatial and temporal periodicities of the optical mesh lattice, there are an infinite number of replicas of the defect straight lines. However, as discussed in the paper, for a sufficiently slowly varying defect, the corresponding Fourier-space defect potential is highly peaked around certain momenta for each line and becomes small elsewhere; this relative importance is schematically represented by an increased size of the blue points, with the largest weight being concentrated around $k=0$ and equivalent points in the extended Floquet-Bloch scheme. (c, d) The dispersion as in (a) and (b) respectively, but only with the fundamental $r=r^{\prime }=0$ defect line (now unadjusted for the relative defect weight), in order to highlight the first intersections between these lines and the Bogoliubov bands as labeled by numbers and letters, respectively. (e, f) The logarithm of the absolute value of the differential intensity spectra [Eq. (35)], as calculated at the late time step $m=900$ for a defect with spatial width ${\sigma }_n=0.25$ (black line) and ${\sigma }_n=0.75$ (red line) moving through a uniform light field with $\mathrm{\Gamma }{I}_0=0.25$ at the speed values (e) $s=1.5v_{\mathrm{sound}}$ and (f) $s=0.5v_{\mathrm{sound}}$. The defect parameters are ${\sigma }_m=100, m_d=500$ and ${\phi }_0=0.01$. Labeled arrows indicate the momenta of the lowest intersections shown in (c) and (d) back-folded into the range $-\pi /2\le k<\pi /2$.

Figure 17

Critical velocity as calculated from the Bogoliubov dispersion within a region $-k_c\le k\le k_c$. For $k_c=2\pi$ or larger, $v_c=0$ for all stable nonlinearities, indicating that the system is not superfluid. However, for lower values of $k_c$, a nonvanishing $v_c\ne 0$ is found in the presence of a finite nonlinearity, indicating that “superfluid-like” behavior may be observed for sufficiently wide defects.

Figure 18

(a) Spatiotemporal colorplot of the numerical differential intensity [Eq. (33)] for a temporally smoothly varying, spatially narrow, and very weak defect (${\sigma }_m=200, {\sigma }_n=1, {\phi }_0=0.00001, m_d=500, n_d=0$), moving at speeds $s=\{1,1.5,2\}{v}_{\mathrm{sound}}$ as indicated by the black parallelogram across a uniform light field, with nonlinearity $\mathrm{\Gamma }{I}_0=0.1$. (b) The total emission intensity [Eq. (49)] plotted as function of the defect velocity and nonlinearity for $m_{\max }=1000$ and $j=101$. As a guide to the eye, we have marked the speed of sound (black dashed line); lower and upper critical velocities for the first Bogoliubov band (red dotted line and blue dashed line, respectively); the lower critical velocity for the second Bogoliubov band (yellow dashed line) within $-\pi \le k,\theta <\pi$. (c) The fraction of states [Eq. (50)] in the Bogoliubov bands which approximately satisfy the restricted Landau criterion as calculated with $N_T=15000$ and $\delta \theta =0.001$. Other marked lines are the same as in panel (b).

Figure 19

(a) Spatiotemporal 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}_{\mathrm{sound}}$ across a uniform light field, with nonlinearity $\mathrm{\Gamma }{I}_0=0.25$, and with all other parameters as in Fig. 8. The smooth defect profile (indicated by the black parallelogram) predominantly excites sound waves, with the strongest emission of radiation for $s=v_{\mathrm{sound}}$. (b) The total emission intensity [Eq. (49)], plotted as function of the defect velocity and nonlinearity for $m_{\max }=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 [cf. Eq. (48)].

Figure 20

(a), (b) Spatiotemporal colorplot of the numerical differential intensity [Eq. (33)] for, respectively, a weak defect with amplitude ${\phi }_0=0.01$, and a strong defect with amplitude ${\phi }_0=0.2$. The defect (indicated by the black parallelogram) moves at a speed $s=v_{\mathrm{sound}}$ in a uniform light field, with nonlinearity $\mathrm{\Gamma }{I}_0=0.45$, and with all other parameters as in Fig. 8. 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 $m_{\max }=1000$ and $j=101$. The instability appears at large nonlinearities, spreading from the regime, $\mathrm{\Gamma }{I}_0>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 $m_{\max }$.

Figure 21

Total emission intensity (49) for a narrow and rapidly varying defect, moved across a Gaussian light field with initial width ${\sigma }_G=10$ (top panel) and ${\sigma }_G=5$ (bottom panel). The defect parameters are ${\sigma }_n=1, {\sigma }_m=10, m_d=20$ and ${\phi }_0=0.01$. The total emission intensity (49) has been evaluated with $m_{\max }=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.