Towards the experimental observation of turbulent regimes and the associated energy cascades with paraxial fluids of light

The propagation of light in nonlinear optical media has been widely used as a tabletop platform for emulating quantum-like phenomena due to their similar theoretical description to quantum fluids. These fluids of light are often used to study two-dimensional phenomena involving superfluid-like flows, yet turbulent regimes still remain underexplored. In this work, we study the possibility of creating two-dimensional turbulent phenomena and probing their signatures in the kinetic energy spectrum. To that end, we emulate and disturb a fluid of light with an all-optical defect using the propagation of two beams in a photorefractive crystal. Our experimental results show that the superfluid regime of the fluid of light breaks down at a critical velocity at which the defect starts to exert a drag force on the fluid, in accordance with the theoretical and numerical predictions. Furthermore, in this dissipative regime, nonlinear perturbations are excited on the fluid that can decay into vortex structures and thus precede a turbulent state. Using the off-axis digital holography method, we reconstructed the complex description of the output fluids and calculated the incompressible component of the kinetic energy. With these states, we observed the expected power law that characterizes the generated turbulent vortex dipole structures. The findings enclosed in this manuscript align with the theoretical predictions for the vortex structures of two-dimensional quantum fluids and thus may pave the way to the observation of other distinct hallmarks of turbulent phenomena, such as distinct turbulent regimes and their associated power laws and energy cascades.


Introduction
The concept of fluids of light lays on a mathematical analogy between the model describing a quantum fluid and that for a laser propagating on or exciting a nonlinear medium [1,2]. The major interest in these analogue systems resides in the fact that they constitute versatile testbeds to observe experimentally quantum-like phenomena and behaviors that are otherwise untestable [3][4][5][6]. Indeed, featuring tunable properties and parameters, they allow to probe a wide range of quantum-like dynamics with inexpensive tabletop optical setups, establishing an interesting landscape of opportunities that have been extensively explored in recent years [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. While most of the experimental work has so far focused on the observation of superfluidic-like signatures in perturbative regimes [9,13,16], there has been an emergent interest in realizing and probing turbulent regimes and dynamics with these analogue systems [12,[18][19][20][21].
To this day, turbulence remains one of the open challenges of fluid and quantum-fluid dynamics. The major obstacle is that turbulence involves the interplay between nonlinear and non-equilibrium dynamics and a large number of degrees of freedom, which rapidly turns conventional theoretical approaches impracticable [22]. Yet, while a complete analytical description is unavailable, it is widely accepted that turbulence reveals statistical properties and signatures shared amongst diverse physical systems and limits, such as vortex excitation and decay [20,[23][24][25][26], and energy cascades with distinct power-laws [23-25, 27, 28]. The hypothesis of the existence of a such class of universal behavior fostered many works in recent years that aim at identifying and probing distinct regimes of turbulence, from weak and wave turbulence [25,[29][30][31][32], to strong-vortex-turbulence regimes [20,[23][24][25][26][32][33][34][35]. Still, the accessibility and degree of control of the most common experimental systems-usually classical fluids [22,36] or quantum fluids like Bose-Einstein condensates [20,21,27,35]-are often limiting factors. Within this context, fluids of light appear as an interesting path for exploration with novel opportunities in the field [12,18,19,37,38].
Fluids of light are usually deployed in two distinct configurations, the cavity configuration [1] or the propagation or paraxial geometry [2]. Focusing on the latter, the fluid of light analogy is constructed over the similarity between the mathematical model for the dynamics of the laser along the propagation direction in nonlinear media with that for the time evolution of a Bose-Einstein condensate under a mean-field approximation. In comparison with the cavity geometry, where a polariton fluid is driven by an incident laser, the behavior loses its true quantum nature but still keeps some signatures of a quantum-like system. The major advantage is that it enables analogue systems to work at room temperature without the limitation of dissipation. As referred to before, the study of superfluidic-like signatures is the most relevant achievement of the field, with the most significant highlights being the observation of dissipation-less flows [2], drag-force suppression [9], persistent currents [7] and the creation of coherent turbulent structures [12]. Still, most of the research has so far focused on the weak perturbative regime and described within the Bogoliubov theory for elementary excitations, leaving stronger nonlinear dynamics a subject that is still underexplored, such as the case of turbulence.
Within this context, this article presents an experimental approach to the use of fluids of light in a paraxial configuration to probe regimes of turbulence, aiming to establish a tunable platform for exploring the connection and crossover between weak and strong turbulent dynamics. For this purpose, we focus on an experimental setup in which two lasers at distinct wavelengths propagate inside a photorefractive crystal, creating an analogue two-dimensional fluid with light that is perturbed by a light-induced obstacle. We show that in specific conditions that match those predicted by the theoretical model, the superfluid-like regime is broken and a turbulence regime can be excited. Using a holographic technique, we reconstruct the output complex states and observe signatures of turbulent dynamic regimes, such as vortex excitation and the characteristic power law in the incompressible energy spectrum. Also, to further confirm the results, we demonstrate that the observations match those obtained using numerical simulations, suggesting that these analogue quantum simulation platforms can be an interesting tool in fluid dynamics for the controllable experimental exploration of a wide range of turbulence regimes.

Physical model for the two-dimensional paraxial fluid of light
The central goal of this work is to explore and validate the use of paraxial fluids of light for probing turbulent regimes. For this purpose, we focus on the analysis of an experimental implementation, similar to the one implemented by Eloy et al [12], that explores the dynamics of the interplay of two laser beams at distinct wavelengths, propagating inside the photorefractive crystal as schematically represented in figure 1. In this setup, the first laser E f = E f (r ⊥ ) exp [i(k f z)]e p at λ f = 532 nm acts as the two-dimensional light fluid. A second laser, with is set at a distinct wavelength λ o = 632 nm to avoid coherent effects and induces an all-optical obstacle that interacts with the fluid through cross-phase modulation effects.
The polarization of both beams, given by the vector e p , is aligned with the c−axis of the photorefractive crystal. Then, neglecting the anisotropic response of the crystal and assuming the paraxial approximation, the propagation dynamics of the envelope functions of both light beams are described by the system of coupled equations, where n e is the extraordinary refractive index, k f = 2π/λ f and k o = 2π/λ o are the wavenumbers of the beams, α denotes the absorption, approximated as equal for both beams, and ∆n max = 1/2n 3 e r 33 E ext . For a photorefractive crystal, the saturable nonlinear refractive index variation in the steady-state [12,39] can be given by where are the local intensities for both beams and I sat is a saturation intensity for the crystal. To further simplify the analysis, equations (1) and (2) can be rewritten in dimensionless units by introducing the new spatial variables x ′ = k f √ ∆n max n e x and z ′ = k f ∆n max z, and normalizing the beams envelopes to the intensity of the beam emulating the fluid, Dropping the primes for simplicity, the model can now be expressed as where I ′ sat = I sat /I f and α ′ = α/k f ∆n max . The constant γ = λ f /λ o was introduced to take into account the difference in the wavelengths.
To introduce the fluid of light analogy we can consider the equation (4) that describes the light fluid and apply the Madelung transformation where ϕ is the light fluid phase profile. It is easy to show that we can obtain the Navier-Stokes equations where the last term on the right side of equation (7) is the Bohm potential, known to be related to quantum-like effects. Thus a clear analogy between light and a fluid is established, with the local intensity being associated with the fluid density, while the fluid velocity is related to the phase as v = ∇ ⊥ ϕ. The velocity of the fluid is thus controlled experimentally by selecting the incidence angle of the E f beam in relation to the propagation z-axis. The Bogoliubov dispersion relation is obtained by linearizing equations (6) and (7), which, in the initial units, gives From this, follows the analogue sound velocity of the fluid and the healing length of typical length scales in the transversal plane where ∆n ≈ ∆n max I f /I sat is the refractive index variation away from the optical defect in the Kerr regime.

Breakdown of the superfluid regime and the development of turbulent states
In line with the Landau condition, the existence of a non-null sound velocity suggests that the light fluid can behave as a superfluid for velocities below c s . In this regime, the intensity modulations due to diffraction are suppressed. At speeds greater than the speed of sound, the superfluid regime is broken and the fluid exerts a force on the defect. This force is given by and can be estimated by considering the fluid of light moving with velocity v, and a static obstacle with a [40][41][42]. Furthermore, we also neglect the absorption of the medium and consider that the nonlinear response is in the Kerr regime (I f ≪ I sat ), such that ∆n total ≈ −I f /I sat . Above the critical velocity, the fluid of light may be written as ψ f = ψ 0,f + δψ f , where ψ 0,f corresponds to the uniform component of the fluid and the δψ f represents the perturbations induced by the optical defect on the fluid. Considering only first-order perturbations, we have that where the first term is null by symmetry considerations, and only the perturbations contribute to the force. The different components in equation (12) can be written in terms of its Fourier transforms which, after integrating in r ⊥ and k ′ , leads to Taken into account first order perturbations in equation (4), yields [40] The force can be divided into two components, one parallel to the velocity and another one perpendicular, F || and F ⊥ , respectively. The integration is then carried out by considering the poles of equation (15) and gives [40] where F 0 = π 2 σ 4 V 2 0 ψ 2 0,f , I n (x) are the modified Bessel functions of the first kind, and k c = 2 v 2 − c 2 s . This expression is true for obstacles with size σ ∼ ξ, where the breakdown of the superfluid regime occurs for velocities above the speed of sound. However, in the case of large defects, σ ≫ ξ, this breakdown can occur for velocities below the speed of sound due to the local variation of the fluid density [43]. Thus, assuming a critical velocity v c < c s , we modify equation (16) Above v c , the fluid starts to exert a force on the obstacle along the direction of the velocity direction, which beyond the linear perturbative regime can also cause hydrodynamic instabilities to form around the defect. For velocities just above the critical velocity, these instabilities can result in the emission of pairs of vortices or other characteristic spatial modulations, which can be interpreted as a dissipation mechanism of the fluid. This regime precedes a turbulent state, for which vortex structures are one of the hallmarks.
While turbulence involves a complex interplay of chaotic phenomenology that is hard to track and model in real space, it is still possible to investigate this regime in terms of characteristic statistical signatures of its energy spectra. A common way to observe these signatures is to look at the density-weighted velocity field, divided into its compressible and incompressible components, that satisfy ∇ · u I (r ⊥ ) = 0 and ∇ × u C (r ⊥ ) = 0. As the signatures of different turbulent regimes appear in the momentum space [23,25,44], it is convenient to transform these velocity fields into its Fourier transforms [44]. In this space, the spectral kinetic energy density of the fluid is also decomposed into its incompressible and compressible components E k (k) = E I k (k) + E C k (k), obtained with the expressions This numerical procedure allows to compute the incompressible component of the kinetic energy of the fluid of light at the output of the crystal for different fluid velocities, and thus probe the presence of the characteristic energy power laws. For velocities just above the critical velocity, one of the most typical structures generated in the non-perturbative regime is the vortex dipole, nucleated past a defect. For this situation, and considering a two-dimensional fluid, the limits of the incompressible component of the kinetic energy for low and high k values can be obtained analytically. According to Bradley and Anderson work in [45], these two regimes are associated with distinct power laws, respectively where d is the separation between the two vortices, and The k 1 and k −3 power laws are associated with the far field velocity distribution and the structure of the vortices, respectively, and the transition between the two occurs around kξ ∼ 1 [45]. Since it is known that these two power laws are associated with non-perturbative regimes that usually precede turbulent dynamics, the goal of this work is to validate experimentally these signatures in our two-dimensional fluid of light.

Experimental setup
In this work, our nonlinear medium is a Strontium-barium niobate-SBN:61 crystal with dimensions 5 × 5 × 20 mm 3 and doped with cerium at 0.002%. A simplified scheme of the experimental setup is depicted in figure 1, using two laser beams in the experience. To create the fluid of light we use a collimated Gaussian beam with approximately 615 µm 1/e 2 radius at λ = 532 nm and I f ≈ 8.4 mW cm −2 . The fluid velocity is controlled by applying different grating levels through a spatial light modulator with its image plane at the crystal entrance. As a result, this effectively alters the angle θ at which the beam enters the crystal, which is related to the analogue fluid velocity (v = θ ne k f ). To create the optical obstacle, we use a He-Ne laser at λ = 633 nm, that passes through a spatial light modulator to create a droplet beam [46]. Due to interference between the two Bessel beams utilized to form the droplet, this beam resembles a localized Gaussian potential and lacks the typical intensity modulations around the central maximum presented in a Bessel beam [46]. We take into account these droplet beams to remove the unwanted effects often produced by the Bessel secondary lobes near the defect, which can impact the energy calculations. Typically, the obstacles induced in our experiments have a diameter of D o ∼ 120 µm and P o ∼ 68.5 µW. The output plane is imaged by a 4f-lens system with a magnification of 2.5 × in a CMOS camera. Additionally, by interfering the first laser with a reference arm, we can obtain an interferogram to qualitatively verify the presence of vortices. This configuration also allows using an off-axis holographic technique to reconstruct the complex description of the beam (profile and phase) [47]. This is required for calculating the kinetic energy and characterizing the turbulent states and will be fundamental for the observation of the predicted power laws.
Lastly, the crystal was previous characterized by using the method presented in [48], and we estimate to have a ∆n ∼ −0.3 × 10 −4 with an I sat = 47 mW cm −2 . Furthermore, the absorption was also estimated to be α ≈ 0.116 cm −1 for both beams. These values were obtained for an applied static electric field E 0 = 1400 V cm −1 and with n e = 2.36.

Numerical simulations
To better understand the dynamics of the experimental setup, we used a numerical solver based on General Purpose Graphics Processing Unit computing [49,50] that integrates the complete coupled system described by equations (4) and (5). It should be noted that despite the model neglecting the anisotropic nonlinear response of the crystal, still yields predictions consistent with the experiments. To make the simulations more accurate, we used the fluid of light profile measured at the crystal input ( figure 1(b)) as the initial condition. For the optical obstacle, we simulated the evolution of a droplet beam in free space and used the droplet state that was most similar to the profile at the crystal entry. In figures 2 and 3, we have different comparisons between the numerical and the experimental results, and we see that they are generally in agreement.
Having a numerical solver to simulate the experimental setup is a useful tool to optimize the experimental parameters to emulate specific phenomenology and explore new configurations that are not yet supported by our setup.

Results
To explore the emergence of turbulent states, we evaluated the fluid output for various fluid velocities. First, we measured the drag force on the defect to observe the breakdown of the superfluid regime. Then, by reconstructing the fields and computing the associated incompressible component of the kinetic energy, we characterize the different regimes by identifying their characteristic power laws.

Drag-force
For every fluid velocity studied, we measured the intensity profiles of the fluid and the optical defect at the output of the crystal. With this data, the drag-force was calculated using equation (11) and the results are plotted in figure 2 along with the results for the numerical simulations and the fit to the theoretical curve given by equation (16). An agreement between the analytical model and the obtained data validates our experimental approach and aligns with the phenomenological explanation introduced in section 2.1.
The superfluid regime is observed for small velocities when the force on the obstacle is null. As the velocity increases, the superfluid regime is broken and a force appears on the defect. Note that the velocity at which this breakdown occurs is below the sound velocity, v c < c s , since we are considering an optical defect much larger than the healing length, specifically D obs ∼ 12ξ. After reaching a maximum value, the system enters the supersonic regime, and the drag force decreases. In this situation, the fluid is sufficiently fast to pass through the defect without causing a large interaction. The theoretical curve given by equation (16) agrees very well with the experimental and numerical data. Note that the force was derived for a static Gaussian defect, and both in the experiment and simulations the optical defect and the light fluid evolve interdependently.
Above the critical velocity, the existence of this force indicates that the system enters a dissipative regime where energy is transferred from the optical obstacle to the light fluid. This exchange is also the main mechanism to excite intensity modulations or nucleate quantized vortices. In the next section, we look at the kinetic energy of the different outputs and analyze the existence of power laws associated with the generated turbulent structures.

Spectrum of the incompressible kinetic energy
The output states, i.e. the amplitude and the complex phase, were reconstructed using the off-axis digital holography method [47]. This procedure allows us to determine the incompressible component of the kinetic energy using equation (19), which is plotted in figure 3. The off-axis digital holography technique was also applied to the numerical data to more accurately approximate the simulations with the experimental data. Note that this approach naturally filters out some of the noise and the output velocity, allowing the study to focus on the turbulent structures created by the obstacle. Furthermore, in order to help the qualitative analysis, we also normalized the output states, such that´|E f | 2 dr ⊥ = 1.
For fluid velocities below the critical velocity, v < v c , the fluid is in the superfluid regime and consequently no characteristic power laws appear, as illustrated in figure 3(a). In this regime, the presence of the defect almost has no impact on the flow of the fluid. Indeed, the intensity modulations visible in figures 3(d) (experimental) and 3(j) (simulation) result from the fact that the state at the input is not exactly the stationary state corresponding to the superfluid flow around the defect, thus producing shock waves. While the experiments are limited to the length of the crystal, the simulations show that these modulations disperse after some propagation distance as the initial state eventually decays to the ideal superfluid state.
For fluid velocities above the critical velocity, v > v c , the flow experiences the nucleation of vortex pairs of opposite vorticity in the wake of the defect, constituting a vortex dipole, as shown in figures 3(f) (experimental) and 3(g) (simulation). In this situation, both experiment and simulations exhibit the signature power law k −3 for kξ > 1 in the spectrum, as depicted in figure 3(b). However, at lower wavelengths, kξ < 1, the expected k 1 power law cannot be identified. An explanation for this absence may reside in the fact that this power law is essentially determined by the velocity field generated by the vortex dipole in the far field. As such, it is sensitive to the relative position of the two vortices, as well as to noise and additional perturbations from the initial shock waves. Furthermore, the integration area for calculating the energy is also limited, thus reducing the sampling of large length scales associated with this power law. As a result of the combination of these factors, the k 1 power law in the spectrum is probably underrepresented, and thus only the transition to the k −3 regime is observed [45].
Finally, for velocities much higher than the critical velocity, v ≫ v c , the distance between the two vortices reduces considerably, as shown in figures 3(h) (experimental) and 3(n) (simulation), suggesting that they may asymptotically merge to form a snake-like instability or a dark soliton [12], or a caviton [51,52]. The merging of the two vortices is difficult to confirm experimentally due to the limitations in the reconstruction method used to recover the output state when the vortices are very close. In this configuration, the power law k −3 is still identifiable in figure 3(c) but over a small range of values of k, while again there is no evidence of the power law k 1 .
The limited size of the crystal allows only to excite one vortex dipole and observe the associated power law. In a perfect situation, where the crystal could extend into infinity along the propagation direction, the defect would continue to pump energy into the fluid, resulting in the nucleation of new vortex dipoles. In this case, the light fluid would develop into a more complex turbulent state where other power laws, hallmarks of different energy cascades, could be seen. Additionally, it would be possible to investigate the transition between weak and strong turbulent regimes by examining the relationship between the compressible and incompressible components of the kinetic energy.

Discussion and concluding remarks
The quest for a better understanding of the phenomenology and mechanisms behind turbulence has long pursued the capacity to generate, control, probe and capture the dynamics of controllable physical systems capable of experiencing a multitude of turbulent regimes and structures. Lacking a complete analytical description, the experimental determination of statistical signatures is a common path to characterize turbulent phenomena and distinct regimes. One of the most popular approaches is to look at the energy spectra and seek for the existence of characteristic scaling power laws, that can be compared with predictions of theoretical models for their validation and improvement. Within this context, the subject of fluids of light in paraxial configurations may provide an interesting paradigm and a wide range of opportunities for simulating and investigating turbulent dynamics under tunable conditions in truthful two-dimensional setups.
In this manuscript, we investigated the capacity of using these analogue systems to emulate turbulent regimes and observe the signature power laws that are theoretically predicted to occur. For this purpose, we created experimentally a fluid of light perturbed by an all-optical defect using two optical beams propagating inside a photorefractive medium and interacting through cross-phase modulation. To provide quantitative evidence of the transition between a dissipation-less regime to a dissipation one, we analyzed the drag force exerted by the fluid on the optical defect, identifying the existence of a superfluid regime below a critical velocity, in agreement with theoretical and numerical models. Furthermore, by using the off-axis digital holography, we recovered the output states of the fluid and determined the spectra of the incompressible component of the kinetic energy. The results show the existence of a characteristic power law k −3 towards smaller scales, which to our best knowledge, constitutes its first experimental observation with paraxial fluids of light.
Despite the simplicity of the structures formed in the turbulent flow generated by the defect for v > v c , this work clearly demonstrates that it is possible to observe the power laws in these flowing dynamics using off-axis digital holography. The next natural step is to apply this technique to more complex systems with more vortices, for example by considering multiple potentials [33], the decay of unstable vortex latices [23,34], or by using the interaction between two fluids of light [37,38]. Hopefully, the states of these systems may allow observing additional signature power laws of various turbulent regimes and their distinctive energy cascades, such as the Kolmogorov k −5/3 power law or the enstrophy cascade with a k −4 power law. This is a research path that we intend to explore in future work.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.