A Eloy

O Boughdad

M Albert

P.-É Larré

F Mortessagne

M Bellec

C Michel

Experimental observation of turbulent coherent structures in a superfluid of light

We experimentally explore the rich variety of nonlinear coherent structures arising in a turbulent flow of superfluid light past an obstacle in an all-optical configuration. The different hydrodynamic regimes observed are organised in a unique phase diagram involving the velocity of the flow and the diameter of the obstacle. Then, we focus on the vortices nucleated in the wake of the obstacle by investigating their intensity profile and the dependence of the radius of their core on the healing length. Our results pave the way for further investigations on turbulence in photon superfluids and provide versatile experimental tools for simulating quantum transport with nonlinear light.

I. INTRODUCTION

Turbulence appears in ordinary fluids as simple as air or water, but also in astrophysical systems [START_REF] Brandenburg | Astrophysical turbluence modeling[END_REF], plasmas [START_REF] Yoshizawa | Plasma and fluid turbulence: Theory and modelling[END_REF] and superfluids [START_REF] Skrbek | Developed suantum turbulence and its decay[END_REF][START_REF] Kolmakov | Wave turbulence in quantum fluids[END_REF]. The main features shared by the different types of turbulence are the unpredictable and irregular character of the system dynamics far from equilibrium, and the interplay between many different physical length scales through highly nonlinear processes [START_REF] Sreenivasan | Fluid turbulence[END_REF]. Quantum turbulence mainly differs from turbulence in ordinary fluids because, in quantum fluids, vorticity is quantised and constrained to appear at a single characteristic length scale [START_REF] Pitaevskii | Bose-Einstein Condensation and Superfluidity[END_REF]. Examples of quantum fluids displaying quantum turbulence are superfluid helium [START_REF] Allen | Flow of liquid helium ii[END_REF][START_REF] Kapitza | Viscosity of liquid helium below the λ-point[END_REF], atomic Bose-Einstein condensates [START_REF] Pitaevskii | Bose-Einstein Condensation and Superfluidity[END_REF][START_REF] Neely | Observation of vortex dipoles in an oblate bose-einstein condensate[END_REF][START_REF] Bloch | Quantum simulations with ultracold quantum gases[END_REF] and quantum fluids of light [START_REF] Carusotto | Quantum fluids of light[END_REF]. Merging nonlinear optics and quantum hydrodynamics, quantum fluids of light have gained great interest in the past few years. Indeed, in properly engineered experimental optical devices, photons can acquire an effective mass and be in a fully controlled effective interaction. They behave collectively as a quantum fluid, and share remarkable common features with other systems such as superfluidity and quantum turbulence. Quantum fluids of light have been investigated in one, two and three effective spatial dimensions (1, 2 and 3D) in various photonic platforms. Among the latter, the major ones are semiconductor [START_REF] Carusotto | Quantum fluids of light[END_REF][START_REF] Deng | Exciton-polariton bose-einstein condensation[END_REF][START_REF] Amo | Superfluidity of polaritons in semiconductor microcavities[END_REF][START_REF] Amo | Polariton superfluids reveal quantum hydrodynamic solitons[END_REF][START_REF] Nardin | Hydrodynamic nucleation of quantized vortex pairs in a polariton quantum fluid[END_REF][START_REF] Sanvitto | All-optical control of the quantum flow of a polariton condensate[END_REF][START_REF] Pigeon | Hydrodynamic nucleation of vortices and solitons in a resonantly excited polariton superfluid[END_REF][START_REF] Pigeon | Sustained propagation and control of topological excitations in polariton superfluid[END_REF][START_REF] Claude | Taming the snake instabilities in a polariton superfluid[END_REF][START_REF] Lerario | Vortex-stream generation and enhanced propagation in a polariton superfluid[END_REF] and optical [START_REF] Kurtscheid | Thermally condensing photons into a coherently split state of light[END_REF][START_REF] Chiao | Bogoliubov dispersion relation and the possibility of superfluidity for weakly interacting photons in a two-dimensional photon fluid[END_REF][START_REF] Chiao | Bogoliubov dispersion relation for a 'photon fluid': Is this a superfluid?[END_REF] microcavities, and cavityless bulk optical media [START_REF] Frisch | Transition to dissipation in a model of superflow[END_REF][START_REF] Pomeau | Diffraction non linéaire[END_REF][START_REF] Leboeuf | Superfluid motion of light[END_REF][START_REF] Carusotto | Superfluid light in bulk nonlinear media[END_REF][START_REF] Larré | Wave pattern induced by a localized obstacle in the flow of a one-dimensional polariton condensate[END_REF][START_REF] Larré | Optomechanical signature of a frictionless flow of superfluid light[END_REF][START_REF] Larré | Propagation of a quantum fluid of light in a cavityless nonlinear optical medium: General theory and response to quantum quenches[END_REF] where the light propagation axis plays the role of time, as sketched in fig.

1.

In the latter propagating geometry (and in a regime where quantum fluctuations [START_REF] Larré | Propagation of a quantum fluid of light in a cavityless nonlinear optical medium: General theory and response to quantum quenches[END_REF] are negligible), fluids of light have been produced in several nonlinear media ranging from liquid crystals [START_REF] Bortolozzo | Optical wave turbulence and the condensation of light[END_REF][START_REF] Laurie | One-dimensional optical wave turbulence: Experiment and theory[END_REF] to thermal liquids [START_REF] Vocke | Experimental characterization of nonlocal photon fluids[END_REF][START_REF] Vocke | Role of geometry in the superfluid flow of nonlocal photon fluids[END_REF][START_REF] Vocke | Rotating black hole geometries in a twodimensional photon superfluid[END_REF] and atomic vapours [START_REF] Šantić | Nonequilibrium precondensation of classical waves in two dimensions propagating through atomic vapors[END_REF][START_REF] Fontaine | Observation of the bogoliubov dispersion in a fluid of light[END_REF][START_REF] Fontaine | Interferences between bogoliubov excitations in superfluids of light[END_REF], but also photorefractive (PR) crystals [START_REF] Wan | Dispersive superfluid-like shock waves in nonlinear optics[END_REF][START_REF] Wan | Nonlinear wave scattering by small barrier potentials[END_REF][START_REF] Wan | Diffraction from an edge in a self-focusing medium[END_REF][START_REF] Sun | Observation of the kinetic condensation of classical waves[END_REF][START_REF] Situ | Dynamics of the berezinskiikosterlitz-thouless transition in a photon fluid[END_REF] in which a quantitative measurement of the normal/superfluid transition has been recently performed [START_REF] Michel | Superfluid motion and drag-force cancellation in a fluid of light[END_REF]. These optical systems allow for a full control and flexibility in the generation, manipulation and observation of photon fluids in various dimensions.

In addition, the nonlinear interactions and the external potentials experienced by the fluid can be precisely customised. Consequently, they consist in versatile systems, fully designated to investigate light hydrodynamics in different kinds of environments, from the simplest case of a homogeneous medium [START_REF] Fontaine | Observation of the bogoliubov dispersion in a fluid of light[END_REF][START_REF] Fontaine | Interferences between bogoliubov excitations in superfluids of light[END_REF] to more complex landscapes with a single obstacle [START_REF] Vocke | Role of geometry in the superfluid flow of nonlocal photon fluids[END_REF][START_REF] Michel | Superfluid motion and drag-force cancellation in a fluid of light[END_REF] or spatial disorder. One of the simplest way to generate turbulence in these systems is to insert an obstacle in the fluid flow. By accurately tuning its velocity, the fluid enters a turbulent regime, giving birth to nonlinear coherent structures ranging from vortex pairs to dark solitons and snake instabilities [START_REF] Claude | Taming the snake instabilities in a polariton superfluid[END_REF][START_REF] Lerario | Vortex-stream generation and enhanced propagation in a polariton superfluid[END_REF][START_REF] Armaroli | Suppression of transverse instabilities of dark solitons and their dispersive shock waves[END_REF][START_REF] Gallemi | Multidimensional josephson vortices in spin-orbit coupled bose-einstein condensates: snake instability and decay through vortex dipoles[END_REF]. Although optical vortices have been widely studied within the last three decades [START_REF] Coullet | Optical vortices[END_REF][START_REF] Shen | Optical vortices 30 years on: Oam manipulation from topological charge to multiple singularities[END_REF], their generation through the process of turbulence have been subjected to few studies, either numerical [START_REF] Frisch | Transition to dissipation in a model of superflow[END_REF][START_REF] Armaroli | Suppression of transverse instabilities of dark solitons and their dispersive shock waves[END_REF][START_REF] De Rossi | Snake instability of one-dimensional parametric spatial solitons[END_REF][START_REF] Coddington | Experimental studies of equilibrium vortex properties in a bose-condensed gas[END_REF][START_REF] Manjun | Controlling the transverse instability of dark solitons and nucleation of vortices by a potential barrier[END_REF][START_REF] Aioi | Controlled generation and manipulation of vortex dipoles in a bose-einstein condensate[END_REF] or experimental [START_REF] Claude | Taming the snake instabilities in a polariton superfluid[END_REF][START_REF] Lerario | Vortex-stream generation and enhanced propagation in a polariton superfluid[END_REF][START_REF] Mamaev | Propagation of dark stripe beams in nonlinear media: Snake instability and creation of optical vortices[END_REF][START_REF] Raman | Vortex nucleation in a stirred boseeinstein condensate[END_REF][START_REF] Gorza | Observation of the snake instability of a spatially extended temporal bright soliton[END_REF][START_REF] Lagoudakis | Quantized vortices in an excitonpolariston condensate[END_REF][START_REF] Kwon | Periodic shedding of vortex dipoles from a moving penetrable obstacle in a bose-einstein condensate[END_REF]. Related works on the Berezinskii-Kosterlitz-Thouless transition and its underlying vortex dynamics have been recently brought to the photonic realm [START_REF] Situ | Dynamics of the berezinskiikosterlitz-thouless transition in a photon fluid[END_REF][START_REF] Small | Berezinskii-kosterlitz-thouless crossover in a photonic lattice[END_REF][START_REF] Roumpos | Power-law decay of the spatial correlation function in exciton-polariton condensates[END_REF][START_REF] Caputo | Topological order and thermal equilibrium in polariton condensates[END_REF][START_REF] Comaron | Non-equilibrium berezinskiikosterlitz-thouless transition in driven-dissipative condensates[END_REF].

In this letter, we propose a systematic experimental study of the different turbulent coherent structures generated in the wake of a 2D fluid of light passing an obstacle in a propagating geometry. We gather our experimental data in a diagram within regions of apparition of the different nonlinear hydrodynamic structures as a function of the fluid velocity and the obstacle diameter, which is not the standard parameter used in related studies [START_REF] Larré | Wave pattern induced by a localized obstacle in the flow of a one-dimensional polariton condensate[END_REF][START_REF] Larré | Optomechanical signature of a frictionless flow of superfluid light[END_REF][START_REF] Hakim | Nonlinear schrödinger flow past an obstacle in one dimension[END_REF][START_REF] Leboeuf | Bose-einstein beams: Coherent propagation through a guide[END_REF][START_REF] Pavloff | Breakdown of superfluidity of an atom laser past an obstacle[END_REF]. We dedicate the last section to a detailed investigation of the intensity profiles of isolated vortices and of the dependence of their radius on the healing length.

II. FLUID REFORMULATION OF LIGHT PROPAGATION

The propagation in the positive-z direction of a continuous laser beam along a nonlinear PR crystal with a z-invariant optical obstacle is ruled in the scalar, paraxial and monochromatic approximations by the following 2D nonlinear Schrödinger equation for the complex envelope ψ f (r = (x, y), z) of the optical field [START_REF] Boyd | Nonlinear Optics[END_REF]:

i∂ z ψ f = - 1 2n e k f ∇ 2 r ψ f -k f ∆n(|ψ f | 2 )ψ f -k f δn(r)ψ f . (1)
The different quantities involved in this equation are defined in the following lines. The analogy with the genuine Gross-Pitaevskii equation for the wavefunction of a quasi-2D weakly interacting Bose-Einstein condensate [START_REF] Pitaevskii | Bose-Einstein Condensation and Superfluidity[END_REF] is a well known fact [START_REF] Chiao | Bogoliubov dispersion relation and the possibility of superfluidity for weakly interacting photons in a two-dimensional photon fluid[END_REF][START_REF] Chiao | Bogoliubov dispersion relation for a 'photon fluid': Is this a superfluid?[END_REF][START_REF] Frisch | Transition to dissipation in a model of superflow[END_REF][START_REF] Pomeau | Diffraction non linéaire[END_REF]. Importantly, the behaviour of the laser beam in the transverse r plane mimics the flow of such a matter fluid past an obstacle, the propagation coordinate z playing the role of time, thus defining a fluid of light [START_REF] Carusotto | Superfluid light in bulk nonlinear media[END_REF]. In this context, the propagation constant n e k f = 2πn e /λ f of the beam of wavelength λ f in the medium of refractive index n e is the optical analog of the atom mass. On the other hand, ∆n(|ψ f | 2 ) < 0, responsible for the defocusing nonlinear response of the crystal, mimics repulsive photon-photon interactions and acts over the whole transverse extension of the fluid of light. In this quantum hydrodynamic description of the laser beam, the optical intensity

I f = |ψ f | 2 corresponds
to the local density of the fluid of light and the gradient of the optical phase arg(ψ f ) gives the local flow velocity v = ∇ r arg(ψ f )/(n e k f ). In the paraxial limit, a uniform initial velocity is given by the uniform phase of a plane wave such as v θ/n e , with θ being the angle between the beam and the z axis (green beam of the bottom arm in fig. 1(a)). Finally, in the last term of eq. (1), δn(r) is a z-invariant refractive-index depletion which is photoinduced by a second laser beam (red beam in fig. 1(a)). This emulates the presence of a penetrable obstacle for the fluid of light. The height and the size of such a potential barrier depends on the intensity I ob and the diameter d of this secondary laser beam [START_REF] Michel | Superfluid motion and drag-force cancellation in a fluid of light[END_REF]. We also define, away from the obstacle where the fluid of light remains unperturbed, an analog healing length ξ = [n e k 2 f |∆n(I f )|] -1/2 , which corresponds to the smallest length scale for density modulations in the transverse plane, and an analog sound velocity c s = (n e k f ξ) -1 = [|∆n(I f )|/n e ] 1/2 for the modulation waves of the fluid of light. Both analog quantities depend on the fluid density, and are then experimentally controlled via the intensity I f . Panels 1, 2 and 3 in fig. 1(b) correspond to three typical density snapshots of the evolution (along the z axis) of the fluid of light (green spots), flowing from left to right (along the x axis) around the obstacle (red spots). The oscillatory pattern and the density drops illustrated in the third panel are particularly representative. The former, reminiscent of interferences between the incident and reflected components of the green beam onto the red one, corresponds, in the hydrodynamics language, to linear or cnoidal periodic excitations radiated far away from the obstacle [START_REF] Khamis | Nonlinear diffraction of light beams propagating in photorefractive media with embedded reflecting wire[END_REF]. The latter is representative of a regime of turbulent flow characterised by the nucleation of vortex pairs [START_REF] Pomeau | Diffraction non linéaire[END_REF], which one studies in this letter.

III. EXPERIMENTAL IMPLEMENTATION

In the experiments, the nonlinear medium considered is a 5 × 5 × 20 mm 3 strontium barium niobate (SBN:61) photorefractive crystal additionally doped with cerium (0.002%) to enhance its photoconductivity. Because of the saturable nature of the PR effect, we take care of always working in a Kerr-like regime (∆n

(|ψ f | 2 ) ∝ |ψ f | 2 )
where the correspondence between eq. ( 1) and the Gross-

Pitaevskii equation is strict (see Supplementary Material for details).

The fluid of light laser beam is a cw source delivering a collimated Gaussian beam with a 460 µm 1/e 2 radius at λ f = 532 nm. The input angle θ is tuned from 0 to ±23 mrad by means of a spatial light modulator (SLM) conjugated with a digital micro-mirror device (DMD) optical system.

The same conjugated devices are used to create the z-invariant obstacle which consists in a diffraction-free Bessel cw laser beam operating at λ ob = 633 nm. The diameter d of the obstacle (defined as the diameter of the first zero of the rotation-invariant Bessel distribution) can be adjusted from 10 to 90 µm. Both laser beams are linearly-polarised along the extraordinary axis to maximise the PR effect. The photo-induced refractive-index modification is of the order of -10 -4 (see Supplementary Material for details).

For the detection part, an imaging system composed of a ×20 microscope objective (NA = 0.4) and of a sCMOS camera allows to get the spatial distribution of the nearfield intensity of the laser beams at the output of the crystal. The interference between the fluid of light beam and a reference beam, recombined before the camera, allows to obtain an interferogram (the fringes are separated by a 2π spatial phase shift). Using the off-axis holography method [START_REF] Leith | Reconstructed wavefronts and communication theory[END_REF][START_REF] Verrier | Off-axis digital hologram reconstruction: Some practical considerations[END_REF], both the amplitude and the phase of the optical fields can be reconstructed.

To summarise, our experimental apparatus makes it possible to i) prepare the initial condition of the 2D fluid of light, ii) design an obstacle with controllable strength and size, and iii) measure the hydrodynamic parameters of the fluid of light, i.e., both its density and velocity, at the output of the nonlinear crystal.

IV. DIAGRAM OF EXISTENCE OF THE TURBULENT COHERENT STRUCTURES IN THE FLUID

The Mach number, defined as M = v/c s , is a crucial quantity that allows to run through all the different hydrodynamic regimes of the fluid of light. In our optical setup, M typically varies from 0 to 3. As shown in [START_REF] Michel | Superfluid motion and drag-force cancellation in a fluid of light[END_REF], for an obstacle diameter of the order of ξ, a crossover from a normal (i.e., nonsuperfluid but stationary) to a superfluid regime occurs at M ∼ 0.4. This deviation from the Landau criterion [69] is actually not surprising. Indeed, this criterion, which predicts a sharp normal/superfluid transition at M = 1, is intrinsically valid for small obstacles and conservative dynamics. Therefore, it cannot be strictly applied to our setup where the optical defect is not of weak amplitude and where photon absorption is non-negligible (about 1.5 dB/cm). Beyond Landau's theory, the critical velocity for superfluidity is smaller than the speed of sound and a nontrivial function of the obstacle parameters (see, e.g., refs. [START_REF] Larré | Wave pattern induced by a localized obstacle in the flow of a one-dimensional polariton condensate[END_REF][START_REF] Larré | Optomechanical signature of a frictionless flow of superfluid light[END_REF][START_REF] Hakim | Nonlinear schrödinger flow past an obstacle in one dimension[END_REF][START_REF] Leboeuf | Bose-einstein beams: Coherent propagation through a guide[END_REF][START_REF] Pavloff | Breakdown of superfluidity of an atom laser past an obstacle[END_REF][START_REF] Albert | Breakdown of the superfluidity of a matter wave in a random environment[END_REF][START_REF] Engels | Stationary and nonstationary fluid flow of a bose-einstein condensate through a penetrable barrier[END_REF][START_REF] Dries | Dissipative transport of a bose-einstein condensate[END_REF] for related studies in 1D), and the very nature of the nonsuperfluid/superfluid transition is modified by the appearance of an intermediate nonstationary turbulent regime. The latter is characterised by the advent of nonlinear coherent structures such as dark solitons undergoing a snake instability before breaking into pairs of quantised vortices of opposite circulations [START_REF] Brand | Solitonic vortices and the fundamental modes of the "snake instability": Possibility of observation in the gaseous bose-einstein condensate[END_REF]. All these structures obviously depend on photon absorption in the present setup while superfluid features have been shown to be robust against dissipation [START_REF] Wouters | Superfluidity and critical velocities in nonequilibrium bose-einstein condensates[END_REF].

The aim of this section is to measure such a phase diagram by systematically varying the Mach number of the flow, v/c s , and the obstacle diameter normalised to the healing length, d/ξ, which is here the parameter used to vary the state of the obstacle, as explained below. Most theoretical studies build a diagram of existence of the different hydrodynamic coherent structures emerging in the fluid by varying the Mach number, as we do, and the height of the obstacle potential [START_REF] Larré | Wave pattern induced by a localized obstacle in the flow of a one-dimensional polariton condensate[END_REF][START_REF] Larré | Optomechanical signature of a frictionless flow of superfluid light[END_REF][START_REF] Hakim | Nonlinear schrödinger flow past an obstacle in one dimension[END_REF][START_REF] Leboeuf | Bose-einstein beams: Coherent propagation through a guide[END_REF][START_REF] Pavloff | Breakdown of superfluidity of an atom laser past an obstacle[END_REF][START_REF] Albert | Dipole oscillations of a bose-einstein condensate in the presence of defects and disorder[END_REF]. In our experimental apparatus, the latter variation would be performed by changing the depth of the refractiveindex modulation δn(r) induced by the obstacle beam. However, the range of accessible values is rather limited. Moreover, changing the value of the optical index imposes to adapt in a synchronised manner the intensities of both the fluid of light and obstacle beams (see Supplementary Material for details). But changing the intensity of the fluid of light beam would impact the values of ξ and c s , and thus dramatically change the global state of the system. In this study we then focus on varying another intrinsic parameter of the obstacle, namely, its diameter d, which has the advantage of not altering the state of the incident fluid. This configuration has also been studied theoretically [START_REF] Huepe | Scaling laws for vortical nucleation solutions in a model of superflow[END_REF][START_REF] Pinsker | Transitions and excitations in a superfluid stream passing small impurities[END_REF], but has, to our knowledge, never been explored experimentally.

The fluid of light and obstacle beam intensities are respectively fixed to I f = 42 mW/cm 2 and I ob = 4P ob /(πd 2 ) = 510 mW/cm 2 at z = 0, where P ob is the input power of the obstacle beam. I f above provides constant values for ξ and c s while I ob is kept constant by simultaneously varying P ob and d. The refractive index depletion corresponding to the obstacle is then fixed at a value δn = -1.7 × 10 -4 = 0.8 δn max where δn max is the maximum refractive-index variation for the used experimental parameters. To be clear, the experimental parameters directly relevant for the study of the hydrodynamics of the fluid of light thus are i) the input angle θ of the beam which generates it, corresponding to the flow velocity v, and ii) the diameter d of the beam producing the optical defect, which controls the effect of the obstacle on the flow. Figure 2 shows typical images of the fluid of light beam intensity distribution (top row) and the corresponding reconstructed phase (bottom row). First, panel (a) corresponds to the superfluid phase. It presents an absence of long-range radiation from the obstacle [START_REF] Amo | Superfluidity of polaritons in semiconductor microcavities[END_REF][START_REF] Michel | Superfluid motion and drag-force cancellation in a fluid of light[END_REF]. In the other panels, we observe optical vortices with their distinctive behaviour: a deep localised intensity drop associated with a phase singularity. To be more precise, one observes in panel (b) a pair of vortices created at the boundary of the obstacle, downstream. They might not be clearly visible in the intensity image, but appear unambiguously in the phase representation. In panels (c) and (d), the same description can be made, both on the intensity and phase distributions, but respectively for two and three pairs of vortices. It is worth noting that the vortices always appear in the wake of the obstacle and in pairs with opposite circulations v • dr = ±2π/(n e k f ), where we recall that v = ∇ r arg(ψ f )/(n e k f ) is the local velocity of the fluid. They arise to a certain extent as in classical hydrodynamics where Karman vortex streets may appear in the wake of moving objects. This is illustrated by panel (g) which is a zoom of the phase image of panel (d). The region of interest corresponds to the two clear localised intensity drops observed in the top panel (d). The stream plot of the phase derivative, i.e., of the velocity (black lines and arrows), indeed indicates the opposite circulations of the vortices. This behaviour is in good qualitative accordance with theoretical and numerical studies reported in the literature (see [START_REF] Engels | Observing the dance of a vortex-antivortex pair, step by step[END_REF] for a pedagogical viewpoint on the subject). In panel (e), one can still picture vortices in the phase, but also an oscillating structure (indicated by the white arrow in the intensity), resulting on the merging of two pairs of vortices, which is often referred as a snake instability in the literature [START_REF] Mironov | Structure of vortex shedding past potential barriers moving in a bose-einstein condensate[END_REF][START_REF] Griffin | Breaking of josephson junction oscillations and onset of quantum turbulence in bose-einstein condensates[END_REF]. In the last panel (f), this snake instability moves away from the obstacle, and forms likely a dark soliton (indicated by the white arrow). As expected, no phase singularity is observed for both the snake instability and the soliton.

A systematic analysis of such images for various v/c s and d/ξ allows us to build, in the spirit of refs. [START_REF] Huepe | Scaling laws for vortical nucleation solutions in a model of superflow[END_REF][START_REF] Pinsker | Transitions and excitations in a superfluid stream passing small impurities[END_REF], a (d/ξ, v/c s ) diagram displaying the different hydrodynamic behaviours that the fluid of light exhibits when flowing past the obstacle. This diagram is shown in fig. 3. Each region is represented by a color (green color scale) and each point corresponds to the apparition of a new distinct feature. For instance, the dark green points, linked by white segments to guide the eye, correspond to the set of flow and obstacle parameters for which one vortex pair suddenly appears by increasing the flow velocity. The squared letters indicate the location on the diagram of the intensity and phase images of fig. 2. Note however that the precise shape and location of the diagram boundaries depend on the crystal length, which sets the total propagation time of the fluid. Indeed, the turbulent structures observed need a finite time (which depends on the fluid and obstacle parameters) to be created, which means that their number will scale with the ratio of the system size to their respective nucleation time.

As expected, this diagram emphasises the existence of three main regimes of transport, namely, a stationary nonsuperfluid, normal regime at high velocity (white), a nonstationary turbulent regime at intermediate velocity, where nonlinear structures such as vortices are abundant (intermediate green colors), and the superfluid regime which is characterised by the absence of long-range excitations (dark green). In addition, few salient features can be commented. First, we clearly observe that as d/ξ increases, the turbulent regime becomes predominant. This regime tends to vanish for low d/ξ, as reported in our previous experiment [START_REF] Michel | Superfluid motion and drag-force cancellation in a fluid of light[END_REF] on the normal/superfluid transition triggered by a narrow obstacle (d ξ). By interpolation, the merging point of the upper (normal/turbulent) and lower (turbulent/superfluid) separatrices is estimated to be reached at a nonzero value of d/ξ for a Mach number v/c s < 1. This is likely due to photon absorption in the medium (about 1.5 dB/cm), in qualitative agreement with previous theoretical works on low-dimensional dissipative condensates [START_REF] Larré | Wave pattern induced by a localized obstacle in the flow of a one-dimensional polariton condensate[END_REF]. A full theoretical understanding of these phase boundaries is however still lacking and is currently subjected to further investigation. Second, for intermediate values of d/ξ, from 9 to 13, the turbulent phase mainly consists in the generation of vortex pairs and both normal (white) and superfluid (dark green) phases still exist at respectively high and low flow velocities. Then, for larger d/ξ, more complex turbulent behaviours appear, involving more vortex pairs, snake instabilities or solitonic structures. For the value of v/c s investigated in the experiment, we were not able to recover a stationary nonsuperfluid phase for a wide obstacle. Note that the three vertical transitions visible in fig. 3 are actually not as sharp as represented on the diagram. Experimental data are missing in these areas and the analysis of these transition is still under investigation.

The diagram of fig. 3 reveals the rich variety of the turbulent regimes appearing in a 2D fluid of light in a propagating geometry. Although these different behaviours were one by one studied both experimentally and theoretically in ultracold-atom [START_REF] Pitaevskii | Bose-Einstein Condensation and Superfluidity[END_REF] and cavity-polariton [START_REF] Carusotto | Quantum fluids of light[END_REF] physics, this is the first time, to the best of our knowledge, that this rich zoology is experimentally explored in its entirety in the context of 2D fluids of light. Here, we gather all the different behaviours in one unique general diagram involving the two main parameters of the constrained flow: the Mach number and the obstacle diameter. It allows us to predict, for a given Mach number and a given obstacle diameter, the nature of the flow (normal/turbulent/superfluid) and the possible turbulent structures that will be supported by the system. Many aspects are still under investigation. In the following, we will focus on a detailed study of the spatial profiles of isolated optical vortices and their dependence on the healing length.

V. QUANTITATIVE ANALYSIS OF AN ISOLATED VORTEX

As already discussed, one manifestation of quantum turbulence in 2D fluids of light is the generation of optical vortices characterised by a vanishing intensity associated with a phase singularity. The intensity profile of such a coherent structure is phenomenologically expected to vary over a length scale close to the healing length ξ [START_REF] Coddington | Experimental studies of equilibrium vortex properties in a bose-condensed gas[END_REF][START_REF] Isoshima | Vortex stabilization in boseeinstein condensate of alkali-metal atom gas[END_REF][START_REF] Konotop | Solutions of gross-pitaevskii equations beyond the hydrodynamic approximation: Application to the vortex problem[END_REF]. Following ref. [START_REF] Pethick | Bose-Einstein Condensation in Dilute Gases[END_REF], the intensity profile of one of our vortices as a function of the radial coordinate r from the core is approximately given by

I v (r) = I 0 (r/ξ) 2 1 + (r/ξ) 2 , ( 2
)
where I 0 = I(r ξ) is the intensity of the fluid of light far from the vortex core. Defining the vortex radius r v at I v (r v )/I 0 = s, where s is an arbitrary parameter smaller than 1, we readily get

r v (s) = s 1 -s ξ. (3)
Equation ( 3) predicts a linear dependence of the vortex core radius on the healing length, and the proportionality coefficient directly depends on the intensity threshold s where the radius is measured.

To experimentally investigate the vortex size as a function of the healing length in a fluid of light, we track and study them under various conditions. More precisely, the fluid of light intensity is varied from 6 to 16 mW/cm 2 to get different values of ξ, ranging from 15 to 30 µm. We performed different sets of measurements, for different values of the obstacle diameter and strength. As detailed in the Supplementary Material, keeping a constant strength for different values of the fluid of light intensity imposes to vary the obstacle beam intensity as well. Here it ranges from 1.5 to 4 W/cm 2 , which corresponds to refractive-index variations from 1.17 to 2.09×10 -4 thus ranging from 0.7 to 0.8 δn max . The obstacle diameter varies from 60 to 80 µm and the fluid of light velocity is fixed to 2×10 -3 . As seen in the previous section both intensity and phase are measured. To extract the vortex sizes from each images, we perform the following processing. First, the vortex is isolated from its environment (e.g., a neighbouring vortex or any other intensity fluctuation). To do so, the position of the vortex core is tracked from the phase singularity and the intensity image is cropped to 25 µm × 25 µm (see fig. 4(a)). Second, we perform a 2π azimuthal integration of the intensity profile. The profile is normalised to the unperturbed intensity I 0 in the absence of obstacle. The corresponding dimensionless radial profile I(r)/I 0 is plotted in fig. 4 data). Then, the radius of the vortex r v is measured at s = I(r)/I 0 = 0.2 (red dashed horizontal and vertical lines in fig. 4(b)). Above this threshold, the profiles of vortices with the highest radii are drastically affected by the environment. The theoretical prediction is represented by the gray dashed line, calculated from eq. ( 2) with ξ = 16.3 µm, the measured experimental value for the vortex pictured in fig. 4(a). All the measured vortex radii are plotted in fig. 4(c) as a function of ξ. The error bars are the standard deviations of the vortex radii measured on the same image. The filled red circle highlights the measure of the vortex radius depicted in fig. 4(a) and we used the correspond-ing healing length to plot the theoretical prediction in fig. 4(b). For comparison, the prediction, r v = 0.5ξ for s = 0.2 in eq. ( 3), is plotted (solid red line). Interestingly, the experimental data are in rather good agreement with this theoretical expectation. Plotting separately the data corresponding to a given experimental configuration (e.g., fixed obstacle strength and diameter) does not exhibit some specific tendency. All experimental points converge to the same behaviour, meaning that, as discussed above, the vortex size only depends on the main characteristic length scale of the nonlinear system, namely, the healing length. Neither the environment (e.g., the obstacle diameter and strength) nor the initial fluid velocity seem to strongly influence the size of the vortex core.

VI. CONCLUSION

In this letter, we have presented an experimental study of the breakdown of superfluidity and of the various regimes of transport in a 2D fluid of light in a propagating geometry. First, we have experimentally built a phase diagram gathering the different types of flow and coherent nonlinear excitations characterising them, in between the stationary normal and superfluid phases at high and low velocity, respectively. Second, we have proposed a quantitative study of the radii of the vortices nucleated downstream from the obstacle, and showed a good qualitative agreement with theoretical estimates. Our experimental platform gives the opportunity to study the rich nonlinear hydrodynamics of a superfluid of light flowing past a single optical defect. However, there is no major experimental restriction to design more complex environments made of randomly distributed obstacles, opening the way to the study of the transport of superfluid light in disordered landscapes. The basic mechanism of the photorefractive (PR) effect remains in the photogeneration and displacement of mobile charge carriers driven by an external electric field E 0 [START_REF] Yoshizawa | Plasma and fluid turbulence: Theory and modelling[END_REF]. The induced permanent space-charge electric field thus implies a modulation of the refractive index of the crystal, ∆n(I, r) = -0.5n 3 e r 33 E 0

I(r)/I sat 1 + I(r)/I sat , (1)
where n e is the optical refractive index and r 33 the electro-optic coefficient of the material along the extraordinary axis, I(r) is the intensity of the optical beam in the transverse plane, and I sat is the saturation intensity which can be adjusted with a white light illumination of the crystal, and E 0 is an external electric field applied to the crystal along the c-axis of the crystal. For a crystal of strontium-barium-niobate (SBN:61), r 33 = 235 pm/V and n e = 2.36. The experimental control over E 0 and I sat allows to precisely tune ∆n(I, r) [START_REF] Brandenburg | Astrophysical turbluence modeling[END_REF].

In the main paper, we write the propagation eq. ( 1), formulated to mimic the Gross-Pitaevskii-like equation, with a nonlinear term, k f ∆n(I f ), representing the interactions and an external potential, k f δn(I ob ) acting as a potential barrier. However, there is a more realistic manner to model the nonlinear interactions and the local depletion of the refractive index. We have to consider the propagation of two coupled beams, the one acting as the fluid of light, and the one acting as the obstacle, in a PR crystal considering a saturable isotropic nonlinearity. The propagation equations thus read:

i∂ z ψ f = - 1 2k f n e ∇ 2 r ψ f -k f ∆n(I f + I ob )ψ f ( 2
)
i∂ z ψ ob = - 1 2k ob n e ∇ 2 r ψ ob -k ob ∆n(I f + I ob )ψ ob (3)
where ψ f and ψ ob are the slowly varying envelopes of the optical fields for the fluid and the obstacle, respectively, and considering the bulk refractive index does not change with the wavelength. The coupling thus comes from the nonlinear term. Indeed, the total nonlinear refractive index of the medium varies with the two intensities I f,ob = |ψ f,ob | 2 as ∆n tot = ∆n(I f + I ob ) = -0.5n with n e = 2.36 and r 33 = 235 pm/V for a SBN:61 crystal with a linear polarisation along the polar axis (c-axis) of Figure 1(a) illustrates the input Gaussian beam, corresponding to the initial intensity profile for the fluid of light. The latter propagates in a medium in which we impose a local index variation ∆n(I ob ) = δn = -1.5×10 -4 . Figure 1(b) represents the total index variation, ∆n tot = δn + ∆n(I f ), reflecting the fact that the fluid of light also contributes, nonlinearly, to the refractive index variation. As a consequence, the maximum variation of the total index of refraction is greater than δn, but in the example, this does not influence the fluid of light. However, a problem appears when the intensity of the obstacle beam makes the nonlinear refractive index variation saturate, which is commonly the case in this kind of experiments. Indeed, in this case, the nonlinear contribution imposed arXiv:2012.13280v1 [physics.optics] 24 Dec 2020 by the fluid itself induces an effective refractive index, ∆n eff , lower than the maximum refractive index, ∆n max , fixed by the obstacle. In order to evaluate the effective refractive index variation felt by the fluid of light, we develop an extension of the usual formula for the nonlinear index variation in a photorefractive crystal of SBN:61 [START_REF] Yoshizawa | Plasma and fluid turbulence: Theory and modelling[END_REF]. We consider that the effective ∆n eff experienced by the fluid can be written as ∆n eff = ∆n(I f + I ob ) -∆n(I f ).

(

) 5
When assuming that the obstacle beam makes the nonlinear index saturate: ∆n(I ob ) = ∆n max , the previous equation reads

∆n eff = ∆n max -∆n(I f ) = ∆n max -∆n max Ĩf 1 + Ĩf = ∆n max 1 1 + Ĩf . ( 6
)
It is then obvious that the effective refractive index seen by the fluid of light decreases as its intensity increases.

If we take the complete expression of eq. ( 5), then we need to plot a map of ∆n eff depending on I f as well as on I ob . To do so, one can develop the expression of ∆n eff such as

∆n eff = ∆n max   I ob /I sat 1 + 2I f Isat + I ob Isat + (I f +I ob )I f I 2 sat   , (7)
whose absolute value is plotted as a parametric plot in fig. 2. In this figure, we represent a contour plot whose different colors correspond to different values of |∆n eff |, with I ob on the vertical axis and I f on the horizontal axis. The white lines correspond to iso-index lines, and the orange dashed horizontal line corresponds to a typical I sat of 400 mW/cm 2 . Inspecting the limits of this expression, we can verify that we have the right intuition:

• if I ob I sat I f , then ∆n eff → ∆n max , which is expected, as in this case, the obstacle beam completely makes the nonlinear index saturate. This corresponds to the upper left corner of the figure.

• if I ob I f and I f I sat , then ∆n eff → 0. This is also expected as in this condition, the effect of the obstacle is overshadowed by the fluid, and if the fluid saturates the medium, the system goes back to linear. This corresponds to the lower right corner of the figure.

In the intermediate cases, it becomes obvious that the fluid beam almost never sees the maximum value of the nonlinear index variation. However, one can manage to optimise this value, and, more important, to keep it constant while varying I ob and I f , following the white lines of fig. 2. This is what is done in the experiments presented in the main paper.

II. BOGOLIUBOV RELATION DISPERSION FOR A SATURABLE NONLINEARITY

In this section, we adapt Bogoliubov's theory of linearised fluctuations to get the dispersion relation of the small modulations of the fluid of light in the nonlinear photorefractive crystal, in the simplest configuration where the fluid is uniform and at rest in the transverse plane. To do so, we start from the following nonlinear Schrödinger equation for the slowly varying wavefunction of the fluid: 1, the nonlinear contribution is negligible, and the sound velocity, corresponding to the slope of the linear part of the dispersion relation at low k ⊥ tends to 0. On the other hand, when Ĩ0 1, the saturation manifests in the fact that the shape of the dispersion relation does not change anymore when increasing Ĩ0 . However, the typical linear and parabolic trends of the dispersion for respectively k ⊥ 1/ξ and k ⊥ 1/ξ remain. Nevertheless, it is important in the experiment that the hydrodynamical parameters (i.e., c s , ξ) still depend on the fluid density, so we take care of never exceeding Ĩ0 = 1.

i∂ z ψ = -

FIG. 1 .

 1 FIG. 1. (a) Sketch of the experimental setup. The narrow red beam creates a z-invariant localised optical defect which acts as a potential barrier in the transverse plane (x, y). The large green beam creates a 2D fluid of light evolving along the z-axis. The input angle θ of the green beam with respect to the red beam is tuned via a spatial light modulator (SLM) conjugated with a digital micro-mirror device (DMD). It controls the transverse velocity of the fluid of light (indicated by the green arrow). Both are propagating through a biased nonlinear SBN crystal and are imaged via a microscope objective on a sCMOS camera. Beam splitters (BS) are used to get a reference beam to build the interferogram and recover both the amplitude and the phase of the fluid of light. (b) Three-step sketch of the evolution of the 2D fluid of light (green) past the obstacle (red). The propagation coordinate z plays the role of time. The third step shows the formation of a typical interference pattern upstream and a vortex-pair generation downstream from the obstacle.

FIG. 2 .

 2 FIG. 2. (a-f) Images of the intensity (top row) and of the phase (bottom row) of the 2D fluid of light at the output of the nonlinear crystal for various values of the Mach number v/cs and of the obstacle diameter d. The fluid of light flows from left to right and the central position of the diameter-varying obstacle is indicated by the red crosses. Each image is 245 µm×245 µm. The intensity images are normalised to their maximum value. (a) Superfluid phase, d = 50 µm. (b-d) Regime of vortex-pair generation. d = 55, 65, 80 µm for (b), (c) and (d), respectively. (e) Regime of snake instability, indicated by the white arrow. d = 80 µm. (f) Solitonic regime. A dark soliton is indicated by the white arrow. d = 90 µm. No phase singularities are associated to the structures observed in (e) and (f). (g) Zoom of the phase image (d, bottom panel). The region of interest corresponds to the two clear localised intensity drops observed in the intensity image (d, top panel). The phase derivative (i.e., the fluid of light velocity field, see the text) stream plot (black lines and arrows) indicates that the vortices of the pair have opposite circulations.

FIG. 3 .

 3 FIG. 3. Diagram of the observed phases of the fluid of light, ranging from normal (white) to superfluid (dark green). The phases in between (green color scale) correspond to the observation of turbulent coherent structures: vortex pairs, snake instabilities and solitons. Circles are experimental data delimiting the different phases. The squared letters in each phase correspond to the images in fig. 2.

FIG. 4 .

 4 FIG. 4. (a) Intensity distribution of an isolated vortex normalised to the maximum value of the image. The white bar corresponds to 10 µm. (b) Intensity radial profile I(r) obtained after an azimuthal integration (see (a)). The profile is normalised by the intensity I0 measured in the absence of any obstacle. The green curve corresponds to the experimental data. The dashed gray curve is the approximate theoretical intensity profile Iv(r) obtained for ξ = 16.3 µm. The vortex radius rv is measured at the threshold s = I(r)/I0 = 0.2. The corresponding value is reported in (c) (filled red circle). (c) Vortex radius versus ξ. The red circles correspond to experimental data. The error bars are the standard deviations of the vortex radii measured on the same image. The red line corresponds to the theoretical law, rv = 0.5ξ obtained for s = 0.2.

  Supplementary materials: Experimental observation of turbulent coherent structures in a superfluid of light A. Eloy, O. Boughdad, M. Albert, P.-É. Larré, F. Mortessagne, M. Bellec, and C. Michel Université Côte d'Azur, CNRS, INPHYNI, France I. PHOTOREFRACTIVE EFFECT AND EFFECTIVE NONLINEAR INDEX VARIATION

FIG. 1 .

 1 FIG. 1. (a) Gaussian initial condition for the fluid beam, (b) effective index variation for the fluid beam, when it modifies itself the refractive index of the material through the photorefractive effect.

FIG. 2 .

 2 FIG. 2. Cartography of the absolute value of the effective nonlinear index variation experienced by the fluid as a function of the fluid intensity on the horizontal axis, and the obstacle intensity on the vertical axis. The saturation intensity is denoted by the orange dashed horizontal line. The white lines correspond to iso-index lines.

1 2kn e ∇ 2 ψ 2 +FIG. 3 . 1 k 2 n 18 )Figure 3

 22312183 FIG. 3. Dispersion relation of the elementary excitations on top of a homogeneous 2D fluid of light at rest. Blue, red and green solid lines, as well as the dotted-black line correspond to the dispersion relation (16) for different values of Ĩ0. The dashed lines correspond to the k ⊥ 1/ξ sonic regime, dominated by the nonlinearity. The slope of the straight lines directly provides the Bogoliubov sound velocity. In the deep saturation regime, the Bogoliubov dispersion relation is independent of the value of Ĩ0, as the superposition of the green and black-dotted curves reveals.

[ 1 ]

 1 O. Boughdad et al. Anisotropic nonlinear refractive index measurement of a photorefractive crystal via spatial selfphase modulation. Opt. Express, 27(21):30360, 2019. [2] C. Denz et al. Transverse-Pattern Formation in Photorefractive Optics. Springer-Verlag, Berlin, 2003.

  3 e r 33 E 0

	Ĩf + Ĩob 1 + Ĩf + Ĩob	(4)

ACKNOWLEDGMENTS

The authors thank S. Nazarenko and S. Rica for fruitful discussions. This work received funding from the European Union Horizon 2020 research and innovation program under Grant Agreement No. 820392 (PhoQuS), the French government, through the UCA JEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01 and the Region Sud.


Rodrigo De Toledo

Bruno Levy

Visualization of Industrial Structures with Implicit GPU Primitives

We present a method to interactively visualize large industrial models by replacing most triangles with implicit GPU primitives: cylinders, cone and torus slices. After a reverse-engineering process that recovers these primitives from triangle meshes, we encode their implicit parameters in a texture that is sent to the GPU. In rendering time, the implicit primitives are visualized seamlessly with other triangles in the scene. The method was tested on two massive industrial models, achieving better performance and image quality while reducing memory use.

Introduction

Large industrial models are mainly composed of multiple sequences of pipes, tubes and other technical networks. Most of these objects are combinations of simple primitives, e.g., cylinders, cones, tori and planes. Each CAD application has its own internal data format and different ways to manipulate the objects. However, the most common way to export data is through a triangular mesh (triangles are currently the lingua franca in all 3D software). Representing simple primitives by tessellated approximations results in an excessive number of triangles for these massive models. Interactive visualization is only possible when applying advanced and complex algorithms such as Far Voxels [START_REF] Gobbetti | Far voxels: a multiresolution framework for interactive rendering of huge complex 3d models on commodity graphics platforms[END_REF].

In this paper we propose a hybrid approach, using GPU ray-casting primitives [START_REF] Toledo | Extending the graphic pipeline with new gpu-accelerated primitives[END_REF], in combination with rasterized triangles, for visualizing industrial models. By substituting GPU primitives (cones, cylinders and torus slices) for most of the original triangles, we can achieve several improvements: quality enhancement (smooth silhouettes, per-pixel depth and shading, and continuity between pipe primitives, see Figure 1); faster rendering speed (CPU/GPU transference reduction and balancing between vertex and fragment pipelines). Memory reduction (we only keep implicit data such as radius, height and position, instead of storing vertices, normals and topology of tessellated models).

Related work

The idea of replacing meshes with other representations for the purpose of interactive visualization has been widely use. Both billboard and nailboard, which is a billboard that includes per-pixel depth information, have been explored in recent GPUs (in the CG tutorial [START_REF] Fernando | The Cg Tutorial: The Definitive Guide to Programmable Real-Time Graphics[END_REF] they are called depth sprites). More complex solutions have been proposed, such as billboard clouds [START_REF] Dcoret | Billboard clouds for extreme model simplification[END_REF] and relief texture mapping [START_REF] Oliveira | Relief texture mapping[END_REF]. In massive model visualization the concept of texture depth meshes was used to obtain interactive rendering [START_REF] Aliaga | Mmr: an interactive massive model rendering system using geometric and imagebased acceleration[END_REF][START_REF] Wilson | Simplifying complex environments using incremental textured depth meshes[END_REF]. Despite being fast, this method presents some visible regions where the mesh stretches into skins to cover missing geometric information. Another issue is storage space (original 482MB PowerPlant uses 10GB of memory).

Recently, a method that has achieved good results for massive models visualization was Far Voxels [START_REF] Gobbetti | Far voxels: a multiresolution framework for interactive rendering of huge complex 3d models on commodity graphics platforms[END_REF]. It is based on the idea of scene partitioning, grouped into a tree by volumetric clusters. Leaf nodes keep the original triangles and inner nodes have a volumetric grid representation. In rendering time, the grid's voxels are splatted, respecting a maximal screen area. The negative points are: very expensive preprocessing; significant memory use (70MB per million vertices); and an overly complex visualization system (including LOD, occlusion culling and out-of-core data management). In an attempt to accelerate the rendering process, most aforementioned algorithms have as side effects reducing quality and/or increasing memory space. As opposed to them, our solution acheives better performance while enhancing quality and reducing memory.

Ray tracing on GPU Since the beginning of programmable GPUs, ray tracing is a target application that explores their parallelism [START_REF] Carr | The ray engine[END_REF][START_REF] Purcell | Ray tracing on programmable graphics hardware[END_REF][START_REF] Wald | Realtime Ray Tracing and its use for Interactive Global Illumination[END_REF]. The z-buffer algorithm is not used to determine visible surfaces because ray tracing already discovers, for each pixel, the closest intersected object given the viewing ray. Purcell et al. [START_REF] Purcell | Ray tracing on programmable graphics hardware[END_REF] broke the ray-tracing algorithm into kernels that access the entire scene data, which is mapped on textures. The idea of encoding the scene in texture was followed by other methods [START_REF] Horn | Interactive k-d tree gpu raytracing[END_REF][START_REF] Popov | Stackless kd-tree traversal for high performance GPU ray tracing[END_REF][START_REF] Günther | Realtime ray tracing on gpu with bvh-based packet traversal[END_REF]. All these methods suffer from GPU memory limitation. They require space to represent geometry and hierarchical structures, such as kD-trees or bounding volumes. In a recent work, Gunther et al. [START_REF] Günther | Realtime ray tracing on gpu with bvh-based packet traversal[END_REF] succeeded in loading the PowerPlant and rendering it at 3 fps.

We also have encoded geometry in texture. However, in contrast with explicitly representing triangles, we record the parameters of implicit surfaces in texture. Thus, we do not have memory space problems (see Section 6.3).

Extended GPU primitives The concept of extended GPU primitives was first introduced by Toledo and Levy [START_REF] Toledo | Extending the graphic pipeline with new gpu-accelerated primitives[END_REF]. They have created a framework to render quadrics on GPU without tessellation. GPU primitives are visualized through a ray-casting algorithm implemented on fragment shaders. The rasterization of a simple proxy triggers the fragment algorithm. To keep GPU primitives compatible with other surfaces, the visibility issue between objects is solved by the z-buffer. Some recent work have concentrated in ray casting cubics and quartics [START_REF] Loop | Real-time gpu rendering of piecewise algebraic surfaces[END_REF][START_REF] De Toledo | Iterative methods for visualization of implicit surfaces on gpu[END_REF].

Topological Reverse Engineering

There are several reverse engineering algorithms for scanned data, composed by a dense set of points laid down on objects surfaces. However, the data we are dealing with present a different situation, since CAD models are generated without any scanning or capturing step. The main differences are: CAD data have sparser samples; the vertices are regularly positioned (spatially and topologically); and it is partially segmented.

Toledo et al. [START_REF] De Toledo | Reverse engineering for industrial-plant cad models[END_REF] has shown that numerical reverse engineering is not the best choice for CAD triangular meshes. Treating sparse data is a weakness for optimization algorithms [START_REF] Thompson | Featurebased reverse engineering of mechanical parts[END_REF][START_REF] Fitzgibbon | High-level CAD model acquisition from range images[END_REF][START_REF] Petitjean | A survey of methods for recovering quadrics in triangle meshes[END_REF]. On the other hand, it is possible to explore the topological information to reconstruct original implicit data. The topological algorithm is exclusively designed to segment tubes composed of cylinders, truncated cones and torus slices ("elbow junctions"). After applying this reverseengineering algorithm on triangular meshes, the result is a set of higher-order primitives (see Figure 2) plus a set of triangles for the unrecovered objects. Two data sets, PowerPlant and P40 Oil Platform, were used for tests (see Table 1).

Storing implicit parameters on GPU memory

In our application, we charge GPU memory with the parameters of the implicit shapes in pre-processing. The goal is to avoid sending this information at each frame to increase frame rate (see results in Section 6). The drawback is that primitive parameters are difficult to modify on the fly. Before the first rendering frame, we include information from all primitives in a floating-point texture, loading it into video memory. Each texel contains four floating-point scalars. We store the following information for each primitive: cone (8 scalars in 2 texels): origin; main axis scaled by height; two radii. cylinder (8 scalars in 2 texels): origin; main axis; radius and height. torus (12 scalars in 3 texels): center; revolution vector; center-to-begin vector; slice angle; two radii.

In texture, the implicit information is grouped by primitive type (Figure 2). We have tested both random and sequential accesses and they do not produce any difference in performance (in NVidia 7900 graphics card) 3 . In rendering time, we draw the primitives by type to avoid too much shader switching.

Recent graphics cards enable the use of floating-point textures that can be read in vertex shaders 4 . In our case, this strategy is better than using vertex buffers because several vertices share the same parameter values. Using a vertex texture makes it possible to represent this level of indirection. This procedure resulted in doubling the speed of our primitives (see Section 6.2).

Rendering implicit primitives on GPU

We have implemented some high-performance GPU primitives aiming at visualization of industrial models. We have developed two specific quadrics extending Toledo's et al. work [START_REF] Toledo | Extending the graphic pipeline with new gpu-accelerated primitives[END_REF]: cones (Subsection 5.1) and cylinders (Subsection 5.2). Based on [START_REF] De Toledo | Iterative methods for visualization of implicit surfaces on gpu[END_REF], we have created a novel torus slice primitive (Subsection 5.3).

Each surface uses a different fragment shader that executes its ray casting. To trigger fragment execution we rasterize the faces of a specific proxy for each primitive type: hexahedron for cones; quadrilateral for cylinders; and an adapted polyhedron for tori. To accelerate ray casting, implicit surfaces are locally described in canonical positions. Fragments that do not intersect the surface are discarded, otherwise the fragment shader computes the z-value (or depth).

The coordinates of proxy vertices encode (u, v) texture coordinate rather than a 3D positional info. For example, we use glRect(-u,-v,u,v) to call the quadrilateral rasterization for one cylinder. The vertex shader reads implicit information from the texture starting in the position (u, v).

Cone sections

The ray casting executed inside the hexahedron uses a local coordinate system. In this system, the base of the cone is zero-centered and it has unit radius and unit total height (distance between the base and the apex). A very simple fragmentshader computes the intersection between a ray and this canonical cone. The fragment shader is also responsible for drawing the cone caps. Note that complete cones (those that include an apex point) are more efficintly rendered by using a pyramid as their bounding-box (see Figure 3).

Cylinders

For cylinders a quadrilateral billboard is the best choice to trigger their fragment shader. It uses only four vertices and it has a very tight projection enclosing the cylinder, which reduces the discarded fragments. The vertex shader computes the vertex position based on the cylinder's dimensions, in such a way that the 4 vertices always form a perfect convex hull for the cylinder's visible body (Figure 4b). The local coordinate system can be deduced directly from the cylinder's main direction and the computed convex hull axes directions. This is a kind of view-dependent coordinate system, where one of the axes is fixed relatively to the world (the cylinder's main axis) and the other axes depend on the viewing direction. Finally, the unique per-vertex information required is the relative position in the convex hull (front/back and left/right). This information is implicitly given by negative/positive values combined with (u, v) texture coordinates.

View-dependent Coordinate System. Generally, in a cylinder's local coordinate system, z coincides with the main axis while x and y are in a plane perpendicular to z, observing the right-hand rule. In our View-dependent Coordinate System, we impose one more restriction: x must be perpendicular to the viewing direction v.

z = (C 1 -C 0 ) | (C 1 -C 0 ) | , x = v × z , y = z × x
The vector v is a normalized vector pointing to the observer (as if it was perpendicular to the sheet of paper for the reader in Figure 4a). Note that y and z are aligned when projected to the screen, although perpendicular in R 3 .

To compute the final position of the vertices P 0 , P 1 , P 2 , P 3 , we must consider the viewer position and perspective distortion. As indicated on Figure 4, we can compute these positions as follows:

P 0 = C 0 + y + x , P 1 = C 0 + y -x P 2 = C 1 + y + x , P 3 = C 1 + y -x , where x = x • r • cos α, x = x • (|O C| -r) • tan α, y = y • r,y = y • r • sin α, and α = arcsin r |O C| .
Thickness control is a special feature developed for the GPU cylinders. In the vertex shader we guarantee that vectors x and x have a minimum size of 1 pixel on the screen. This avoids the dashed line aspect of thin and high cylinder rendering, which may occur when the viewer is far away. Cylinders' caps. Similarly to cones, it is straightforward to compute the intersection between a ray and the canonical cylinder's body. However, to render the cylinder's caps, there are three possible strategies to attach them:

-Extending the billboard in directions P 2 P 0 and P 3 P 1 (Figure 5a).

-Using a partial 3D bounding box, with 2 faces and 6 vertices (Figure 5b).

-Using a second billboard showing the cap turned to the viewer (Figure 5d).

In the first option, note that P 2 P 0 and P 3 P 1 directions are divergent and extending them may include too many fragments that will be discarded at the end. On the other hand, the second and third options have the drawback of increasing the total number of vertices.

In our implementation we have chosen the last option, because it is the one that better reduces the number of discarded fragments. The vertex shader of this extra GPU primitive uses the same parameters stored on the floating-point texture. Notice that only one extra primitive is used to render both caps, because they are never visible together in the same frame. Finally, in our industrial model visualization, most cylinders do not have caps because they are part of a pipe sequence. Thus, the extra billboard is used only for few ones.

Torus slices

In industrial structures, torus shapes appear in junctions, chains and CAD patterns. Actually, in most cases, the objects contain only a slice of a torus. We use a bounding box well adapted to its form rather than a simple hexahedron. This polyhedron was chosen to reduce pixel wasting by better fitting the slice. At the same time, it is possible to surround torus slices from small angles up to 180 degrees without a significant vertex cost.

We have implemented a special vertex shader to automatically locate vertices around the torus slice based on its parameterization. Once more, (u, v) texture coordinates are associated with positive/negative values, determining internal/external and top/base position. The third coordinate represents the relative angle in the interval [0, 1]. The vertex shader fetches torus parameters and computes vertex world and local coordinates.

We use an iterative algorithm for ray casting tori, the Newton-Raphson root finder presented by Toledo et al. [START_REF] De Toledo | Iterative methods for visualization of implicit surfaces on gpu[END_REF]. They have compared several solutions, considering this one the fastest solution for ray casting torus on GPU.

Results

Enhancing Image Quality

There are five image quality improvements obtained with our method.

Smooth silhouette. The tessellation of curved surfaces impedes the rendering of continuously smooth silhouette, since it is restricted by mesh discretization. The GPU primitives are computed by pixel, therefore the silhouettes are always smooth, even after a huge zoom. See Figure 1(a).

Intersections (Per-pixel z-computation). In conventional triangle rasterization, the z-buffer depth of each pixel is the result of interpolation of pervertex depth. On the other hand, the GPU primitives compute per-pixel depth, which is much more accurate. When two or more GPU primitives intersect, the boundary has a correct shape due to this precise visibility decision (Figure 1(b)).

Continuity between pipe primitives. In industrial plant models, the pipes are a long sequence of primitives. With the tessellated solution, there is a risk of cracks appearing between primitives. To avoid them, the tessellation should keep the same resolution and the same alignment for consecutive primitives. In Figure 1(c), a misalignment has caused cracks. This kind of undesirable situation is also a problem in applications using LOD in pipes (see [START_REF] Krus | Adaptive tessellation of connected primitives for interactive walkthroughs in complex industrial virtual environments[END_REF]), whereas with GPU primitives continuity is natural since there is no discretization.

Per-pixel shading. Our primitives compute shading by pixel (Phong shading), enhancing image quality if compared to default Gouraud shading.

Cylinder thickness control. The thickness control adopted in the vertex shader of GPU cylinders avoids dashed rendering when zooming out of a thin and high cylinder. Moreover, when this control is associated with color computation based on the normal average, it significantly reduces the aliasing pattern effect for a group of parallel cylinders (see Figure 6b).

Performance

We have done two sets of performance measures targeting different comparisons. The first set uses the PowerPlant and compares GPU primitives with tessellated data of the original model. The second set uses P40 model to compare the visualization of GPU primitives with different levels of mesh tessellation.

All tests were done on an AMD Athlon 2.41GHz, 2GB memory, with GeForce 7900 GTX512MB. We render the models entirely in frustum view, fulfilling a 1024 × 768 screen. In this scenario, the bottleneck is not in the fragment shader because primitives are not so large on the screen.

Important Note. All the rendering tests were executed without any culling technique, which would surely speed up all the frame rates. We have intentionally presented the results in this way to avoid covering up speed results when comparing GPU primitives and exclusively-rasterized polygons. PowerPlant. PowerPlant sections 1, 15, 19 and 20, used in Table 2, were intentionally chosen because they were 100% recovered by the reverse engineering algorithm. This way, we can directly compare rasterization and GPU primitive techniques. We have created two situations for GPU primitives: isolated (without using the floating-point texture, but passing their implicit information at each frame) and grouped. The latter is twice as fast as the isolated solution. We also have compared two rasterization strategies based on VBO (Vertex Buffer Object): decomposing the model in multiple VBOs and using a single VBO. Decomposing gives better results for large models.

Comparing the fastest solutions for GPU primitives and for triangle rasterization (second and third columns in Table 2), grouped GPU primitives are clearly faster, almost doubling the speed with multiple VBO implementation.

A different situation is presented on the two-last lines of Table 2, where some of the original triangles (10%) were not converted to implicit primitives (e.g., walls). We compare hybrid rendering (GPU primitives + UT, second column) with the best VBO rasterization rendering (third column). The hybrid solution is between 40% and 60% faster than rasterization.

Oil platform -P40. In this set of measures we have compared the visualization of GPU primitives with mesh rasterization with different levels of tessellation.

The tests were done with 4, 8 and 12-sided meshes (see Figure 7). We were not able to load a 16-sided mesh version because of memory restriction.

As can be seen in Figure 7, the GPU primitives achieve a much better image quality than any tessellated solution. The performance of GPU primitives is comparable to the 8-sided solution (see Table 3).

Memory use

In Table 4 we summarize memory-space information before and after recovering implicit surfaces. The topological recovery procedure converts 90% of the original data of industrial models (PowerPlant section 1 is an exception because it is only composed by tubes and pipes resulting in 100% of conversion rate). The memoryspace of recovered data is reduced to at most 2% (98% of reduction). This is a consequence of compact implicit representation used for recovered primitives. If

Model

Triangles Initial space Conversion Primitives space Final space PP Sec. [START_REF] Gobbetti | Far voxels: a multiresolution framework for interactive rendering of huge complex 3d models on commodity graphics platforms[END_REF] we consider the remaining 10% of unrecovered triangles, industrial models such as oil platforms and power plants can be stored in about 15% of their original data size. The size of floating-point textures on recent graphics cards is limited to 4096 × 4096 (16M texels with 4 floating-points), using 256MB. Our biggest example, P40 model, only uses 4% of this space. Note that it is also possible to use multiple textures, which would push further the scene-size limits.

Conclusion

The use of GPU primitives for CAD and industrial models is very promising. The benefits are classified in three categories: image quality (e.g., perfect silhouette and per-pixel depth), memory and rendering efficiency. Grouping primitive information in a GPU texture has proven to be a good strategy for performance purposes without any memory space problem.

As future work, we suggest the application of well-known acceleration techniques that are usually applied on conventional visualization of massive models: frustum culling, occlusion culling, coherent memory cache, and so on. Since GPU primitives update the Z-buffer, special techniques, which are usually applied to triangles, can also be combined with our primitives. Two interesting future work are shadow maps application and occlusion culling based on queries.

Most implementations of GPU ray tracing have adopted triangles as their only primitive. We suggest the use of implicit primitives, which would probably reduce memory use, which is a critical issue for them.

We also recommend the implementation of other surfaces found in industrial plant that were not covered in our work (e.g., sheared cylinders and half spheres).

Fig. 1 .

 1 Fig. 1. An industrial piece in PowerPlant (section 13). Left: Original tessellated data. Right: Rendering with implicit GPU primitives. Image quality enhancement: (a) silhouette roundness; (b) precise intersection; (c) continuity between consecutive primitives.

Fig. 2 .

 2 Fig.2. Recovered primitives: cylinder, cone section and torus slice. Grouping all implicit information in a floating-point texture (each texel has four floats).

Fig. 3 .

 3 Fig. 3. The front faces of the cone's bounding-box are used to trigger the fragment shader responsible for rendering the cone. A billboard is the best solution for a cylinder without caps. For torus slices, we use an adapted bounding polyhedron with 14 vertices. Each vertex has a parameterized coordinate representing base/top, external/internal and partial angle (1 is the complete slice angle, in this case, 90 degrees).

Fig. 4 .

 4 Fig. 4. (a) View-dependent Coordinate System. (b) Billboard in perspective covering the cylinder's body. (c) We reduce the geometric problem to R 2 (O is the observer projected onto the plane defined by the cylinder's main axis). We can compute the vertices position based on cylinder radius r, distance |O C| and unit vectors x and y.

Fig. 5 .

 5 Fig. 5. Three different ways to render the GPU cylinder with caps: (a) extending the billboard used for the body; (b) using two faces and six vertices of a bounding box; or (c) rendering two separated billboards, one for the body and another one for the caps.

Fig. 6 .

 6 Fig. 6. (a) PowerPlant and P40 Oil-Platform models used on performance tests. (b) PowerPlant Section 1, about 60,000 primitives recovered from 3,500,000 triangles. The thickness control prevents undesired visible patterns when cylinders are seen from far.

Fig. 7 .

 7 Fig. 7. Comparing GPU primitives and tessellations for the P40 model.

Table 1 .

 1 Reverse engineering results with PowerPlant and P40 Oil Platform.

	Data set	triangles (∆) cylinders cones tori unrecovered ∆ effectiveness
	PowerPlant	12,742,978	117,863 2,150 82,359	1,251,019	90.18%
	Oil Platform 27,320,034	215,705 40,001 85,707	2,932,177	89.26%

Table 2 .

 2 PowerPlant performance test. In first block, Sections 1, 15, 19 and 20 do not have unrecovered triangle (UT), in other words, triangles were 100% converted into primitives by the reverse engineering. Frames per second (FPS) measured without v-sync. In second block, Section 12 and the entire PowerPlant have some UT.

	Sec.			FPS			Other
		GPU prim.	VBO		information
		isolate group multiple	single #triangles	#UT #prim.
	1	48.5 93.0	51.0	2.9 3,429,528	0 57,938
	15	138.5 268.5 177.5	165.5 1,141,240	0 21,839
	19	62.5 128.0	73.5	3.8 2,650,680	0 42,046
	20	66.5 130.5	66.9	4.1 2,415,976	0 41,872
		group +UT	all only UT #triangles	#UT #prim.
	12	911 719	524	2935	360,872	38,067 6,205
	all	27.8	21	12.9	81 12,742,978 1,251,019 202,372

Table 3 .

 3 Performance results for oil-platform P40. (*) 341,413 primitives: 215,705 cylinders, 40,001 cones and 85,707 tori.

Table 4 .

 4 Memory space comparison.

		3,429,528	119MB	100%	1.9MB	1.9MB
	PPlant	12,742,978	482MB	90.18%	6.5MB	79.3MB
	P40	27,320,034 1033MB	89.26%	10.25MB	121.19MB

Horn et al. [11] indicate that, in their case, coherent access was better than random.

We use the RECT texture target without filtering.