Rayleigh–Taylor instability in nonlinear Schrödinger flow

The Rayleigh–Taylor instability (RTI) is a fundamental fluid instability that occurs when a light fluid is accelerated into a heavier one. While techniques for observing the RTI in classical fluids continue to improve, the instability has not been demonstrated in quantum fluids. Here, we exploit the formal equivalence between condensed matter and coherent nonlinear optics to observe the superfluid-like instability directly in the optical system. For the RTI, an initial refractive index gradient sets the acceleration, while self-induced nonlinear interactions lead to velocity differences and shear. The experimental observations show that density fingering is always accompanied by vortex generation, with perturbation modes following a hybrid dynamics: horizontal modes (along the interface) propagate as an incompressible fluid, but the vertical length scale (mixing length) is set by compressible shock dynamics. The growth rate, obtained analytically, shows that inhibition due to diffraction has the same spectral form as viscosity and diffusion, despite the fact that the system is dispersive rather than dissipative. This gives rigorous support for the observation that turbulence in quantum fluids has the same scaling as turbulence in normal fluids. The results hold for any Schrödinger flow, e.g. superfluids and quantum plasma, and introduce a new class of fluid-inspired instabilities in nonlinear optics.

Stratified flow, which occurs whenever adjacent layers in a fluid have a density or velocity difference, is fundamental to fluid dynamics. For normal fluids, the corresponding instabilities are textbook problems [1]: thermodynamic fluctuations trigger the growth of perturbations, whose spatial scale is set by the resistive forces of viscosity or surface tension. Accelerated interfaces, in particular, plague many classical flows [2,3]. However, the initial metastable state of layered fluids is difficult to set up experimentally, while dissipation makes it difficult to compare measured values with ideal fluid predictions [4]. Superfluids provide a naturally ideal situation, but their transport dynamics are different: the temperature is ideally zero and the coherent/condensed waves experience no viscosity. As a result, instability dynamics for these systems have received far less attention [5][6][7][8][9]. Nevertheless, interface issues will become increasingly problematic for the atomic community as quantum fluids are transported and mixed [10,11].
Linearized Schrödinger theory suggests that viscosity is replaced by diffraction, i.e. the (frictionless) opposition to wavefunction localization [5,6], but nonlinearity is necessary as well, as interactions determine the condensate pressure [12][13][14], coupling between layers [15], etc. Here, we construct an optical system equivalent to a ground-state superfluid [16] that allows the observation of the ideal Rayleigh-Taylor instability (RTI) [17,18] as a function of density (intensity) difference, internal pressure (nonlinearity) and acceleration (refractive index gradient). In parallel, an analytic theory is developed that collapses the parametric studies into a single unified scaling law. The results highlight the role of wave/phase diffusion in quantum transport and give insight into the onset and statistics of superfluid turbulence.

Theoretical background
Perfect quantum fluids, in which all the atoms have condensed into the ground state ψ, may be described by the Gross-Pitaevskii equation [12][13][14]: 3 whereh is Planck's constant divided by 2 , m is the mass of each particle, V is an external potential and γ measures the strength of (mean-field) interactions. This nonlinear Schrödinger equation is an excellent approximation for Bose-Einstein condensates [19] and is often used to describe superfluids when the normal fluid component is not important [20].
Here, we construct an optical system equivalent to a ground-state superfluid by considering a beam propagating in a nonlinear Kerr-like medium [16]. To focus on spatial dynamics only, we consider monochromatic light of a fixed frequency, so that there is no formal time dependence. Paraxial propagation along the axis z is then well described by the nonlinear Schrödinger equation where ψ is the slowly varying amplitude of the optical field, k 0 = λ 0 /n 0 is the wavenumber in a material with base index of refraction n 0 , λ 0 is the wavelength in free space, n(y) is an imposed refractive index profile, and γ measures the nonlinear index change induced by the light intensity (γ < 0 for defocusing). The nonlinear response of the photorefractive crystal used in the experiments below is similar, with the coupling strength γ being controllable by an applied voltage. It is clear that equations (1) and (2) are identical, with propagation along z replacing evolution in t, wavenumber replacing mass and index of refraction acting as a potential. Therefore, the two equations must support the same type of energy flow. A fluid interpretation follows by applying the polar (Madelung) transformation [21] to equation (2), where ρ = |ψ| 2 is the intensity and S is the coherent phase of the wave function. Scaling then gives ∂ρ ∂z +∇ ⊥ (ρu) = 0, This set of equations is a nonlinear eikonal representation, expressing the conservation of intensity (density) ρ and momentum ρu, where u = ∇ ⊥ S is the direction of energy propagation. In optical terms, the pressure P = |γ |ρ 2 /2 arises from a self-defocusing nonlinearity and an effective gravity g = ∂n(y)/∂ y is provided by a refractive index gradient. The last term in equation (4), known as 'quantum pressure' in condensed matter systems, has the highest-order derivatives and regularizes the system. It is an unusual regularization, in the sense that it is a function of density only and yet appears in the velocity equation (in contrast with normal gases and fluids, which have diffusion in the density equation and viscosity in the momentum equation). Without the quantum pressure, equations (3) and (4) represent ideal Eulerian flow, and there would be no limit to the amount of energy that could accumulate in small spatial scales. For example, shock waves would develop infinitely sharp fronts and instabilities would have perturbation growth at arbitrarily short wavelengths.

Dynamics of the Rayleigh-Taylor instability
The basic quantum RTI was considered by Bychkov et al [5], who followed the usual hydrodynamic practice of assuming an initial mechanical equilibrium (obtained by setting the forcing terms in the rhs of equation (4) to zero) and incompressible flow (∇ · u = 0). As diffraction prevents an ideal step initial condition, we consider a smooth, exponential transition region between the two densities/intensities. The resulting growth rate ω for perturbations is given by [5] where α = 1/ρ 0 (dρ 0 /dy) and L is the characteristic width of the interface (e.g. the linear healing/diffraction length). For long-wavelength perturbations, k x /α 1 and k x L 1, the density differential may be approximated by the difference, reducing equation (5) to In this form, the growth rate (5) consists of two parts: the classical Rayleigh-Taylor growth rate [17,18] and a restorative term due to diffraction [5]. Interestingly, the quadratic (k 2 x ) dependence of the resistance is more reminiscent of diffusion [22] than the cubic dependence indicated by equation (4) and is typical of non-dissipative terms, such as curvature from surface tension [1]. This suggests a reason why superfluid turbulence, which can be triggered by the RTI, has the same Kolmogorov scaling law [23,24] as normal fluids (although fluctuations about this mean may be different [25]).
On the other hand, the diffractive form of the result does not include nonlinearity, apparently violating the rule that linear systems are always stable (there is no mode coupling and thus no possibility for wave growth). This 'violation' is an artifact of the (nonlinear) Madelung transformation and simply states that different spatial modes of the initial interface jump propagate (diffract) at different rates. More serious is the presence of dispersive shock waves (DSWs), which always occur across a step of different intensities/densities [16,[26][27][28][29][30][31][32][33]. These prevent any initial mechanical equilibrium (which can be accommodated by translating to a frame in which the interface is moving [34]) and make the assumption of incompressibility quite suspect. Below, we show experimentally that treating α as a dynamic variable self-consistently remedies all of these issues.

Experimental mapping
As described above, the RTI is created by accelerating a fluid into one that is denser. The optical version uses a refractive index gradient to drive together light of different intensities. The experimental setup is shown in figure 1. The high/low intensity input is created by sending 532 nm laser light into a step-index filter and imaging the interface between different attenuation strips. The acceleration is created by passing an ordinarily polarized plane wave through a continuously varying attenuator and imaging it into a photorefractive crystal; due to nonlinearity, the input intensity gradient creates a refractive index gradient in the medium [35,36], so that wherever the light is, it always 'sees' an area of higher index beneath it.
The material of choice is an 8 × 8 × 8 mm SBN:75 (Sr 0.75 Ba 0.25 Nb 2 O 6 ) photorefractive crystal. For this crystal, the nonlinear index change n = −(1/2)n 3 0 r i j E appÎ /(1 +Î ), where n 0 = 2.3 is the base index of refraction, r i j is the appropriate electro-optic coefficient with respect to the applied field E app and the crystalline axes, and the intensityÎ = |ψ| 2 measured relative to a background (dark current) intensity [37]. While this nonlinearity is saturable, the defocusing parameters used here minimize the difference between the photorefractive and Kerr cases [16].
In contrast to the gradient beam, which propagates almost linearly through the crystal, the fluid-like signal wave is extraordinarily polarized and highly nonlinear (r 33 = 1340 pm V −1 versus r 13 = 60 pm V −1 ). Moreover, the internal pressure of the light waves can be controlled by adjusting the self-defocusing bias voltage applied across the crystalline c-axis. The largest nonlinearity corresponds to the largest voltage and intensity, which in our system is n/n 0 = −2 × 10 −4 .
Light exiting the crystal is imaged into a CCD camera. The degree of acceleration is calibrated by measuring the deflection of a low-intensity signal wave. For the experiments with changing voltage, which also changes the strength of the optical induction, the gradient filter is adjusted to maintain a constant acceleration. The velocity (gradient in phase) is identified directly by interfering the output with a reference plane wave.

Experimental results
The dynamics of the interface was tested as a function of all the driving parameters. The basic instability is shown in figure 2. In response to the index gradient and nonlinearity, the interface buckles, creating the interpenetrating fingers characteristic of the RTI. As a test of these factors, calibration runs were performed for the individual parameters in equation (4). In the linear case, without 'gravity', the interface is stable and simply diffracts (figure 2(e)). Adding the index gradient causes a small drift downward (figure 2(f)). In the nonlinear case, without gravity (figure 2(g)), the step difference in pressure leads to a translation of the interface, with intensity oscillations along (but not perpendicular to) the index gradient [16]. These spatially dispersive shock waves occur whenever interactions act in the same direction as dispersion/diffraction and are the dominant mode of transport in dissipationless systems [16,[26][27][28][29][30][31][32][33]. We show below that these shock waves set the characteristic length scale for the mixing layer when instability is present, but by themselves they are relatively stable to transverse perturbations [16]. Only when both the nonlinearity and index gradient are present does an instability form ( figure 2(b)). For the inverse scenario, in which the dense component is beneath the lighter one, the interface is stable ( figure 2(h)). The instability can be characterized by considering the width, spacing and length of the perturbation fingers ( figure 3). In normal fluids, viscosity or surface tension sets the characteristic spatial scale [1]. Initial perturbations are periodic, meaning that the finger width is half their period. Longer evolution times lead to growth and narrowing of the fingers; if the driving force is strong enough, the resulting shear will overcome resistance and create vortices along the edges. In Schrödinger fluids, the lack of viscosity means that nonlinearity is dominant from the very beginning of flow. We observed vortices at the output for every unstable initial condition, accompanied by narrowing of the fingers as the nonlinearity increases [38]. As shown in figure 3(b), the vortices, identified as forked discontinuities in the phase of the coherent wave, outline the edges of the fingers much more distinctly than the intensity pattern ( figure 3(a)). Consequently, all measurements of the characteristic scales are made using this technique.
A particular region of interest is the mixing layer, in which the two (optical) fluids interpenetrate. A direct measurement of the finger lengths ( figure 3(d)) shows that the length follows the nonlinear scaling relation L = a+b √ γρ H /ρ L , where a = L n 3 0 r 33 E app ρ L /2 = 44 µm is a dimensional scaling constant that represents the intrinsic (diffracting or healing) width of the interface and b ∼ 1 (0.98, 0.92, 1.00 as shown in figure 3). This is the same scaling as that found for a step intensity change without acceleration [16,29,30,33], despite the generation and shedding of vortices.
To explore the general features of quantum RTI, we fix the applied voltage across the crystal at −450 V, so that changes in the nonlinearity (interaction strength) are due to the initial intensity difference only. Figure 4(a) shows the characteristic instability period as a function of 9 acceleration and initial intensity jump. As the driving forces become stronger, the perturbations grow with higher energy (higher spatial frequency or smaller characteristic period). At first sight, the spread in data suggests that the variables are independent. However, equation (5) suggests that there is a maximum growth rate, obtained by setting dσ/dk x = 0, corresponding to a dominant mode k max x = (gαL)/α 2 in the dynamics. This dominant mode depends sensitively on the choice of parameters, in particular the strength and characteristic length of initial density difference. Different choices for these parameters are shown in figures 4(b)-(d). Figure 4(b) shows the incompressible limit first derived by Bychkov et al [5]. As discussed after equation (5), the density prefactor is given by the Atwood number, while the decay rate is given by the (constant) diffraction length L = a. There is a clustering of the data, except for the strongest nonlinearity (smallest period), but no true collapse. Figure 4(c) shows the Bychkov result when the measured shock length is used to determine α. In this case, the Atwood prefactor causes the data to separate, with trends no better than the raw data given in figure 4(a). Figure 4(d) shows the hybrid compressible-incompressible limit suggested by the experiments. The density prefactor αL is given by ln(ρ H /ρ L ), which is exact for an exponentially decaying transition region, while the decay rate is taken from the measured shock lengths of figure 3(d). In this case, there is a clear collapse of the data onto a straight line.

Discussion
The analysis above suggests a rather straightforward evolution of the instability: the generation of intensity fingers followed by shear-induced vortex formation along their edges. Initially, noise creates perturbations which have a momentum component in the direction of the index gradient. By Snell's law, these spots of light will bend downwards, into an area of higher index. This will create an even larger angle of incidence, causing more bending, etc. Along the intensity fingers, the relative velocity flow will create vortices due to a superfluid Kelvin-Helmholtz instability [7]. Alternatively, the intensity/density difference between the inside and outside of each finger creates an effective potential well (relative index difference), inducing vortex formation along its boundaries [39].
However, there is another pathway to vortex generation which must be considered: the snake instability of a dark stripe [40][41][42][43]. There is no such stripe initially, as the input phase is constant across the interface, but there are two ways in which it can form: relative slipping between the top and bottom fields and dynamic generation as the leading edge of a spatially dispersive shock wave (DSW). In the former scenario, the two intensities (densities) slip along each other and propagate separately, with the upper one acquiring a faster phase than the lower one. Eventually, there will be a phase difference between them, creating a dark stripe at the interface due to destructive interference. We note, however, that the snake instability works best at high nonlinearity, when both the upper and lower branches have the same intensity [41]. For the conditions here, snaking occurs very slowly. To check this, we removed the index gradient and introduced an initial phase difference between the branches (not shown); even for equal intensities, with the nonlinearity at its highest value, we saw no evidence of kinks.
The phase-slipping scenario is also unlikely because we are in a gravity-dominated regime (by design), which makes motion normal to the interface the most significant driving factor. On the other hand, the interface has the step initial condition of a shock wave, whose oscillating front resembles a train of dark solitons. While the transverse instability of DSWs can follow from the dynamics of their dominant soliton [44,45], simulations of transonic stripes has shown that the evolution of the modulations can be vortex-free [46]. Further, there are many situations, such as tunneling [47], in which the full structure of the DSW must be considered. Indeed, experiments on DSW propagation and collisions performed using the same crystal as the one here have suggested a surprising robustness, if not absolute stability, of shock waves without acceleration [16]. For this reason, and the success of the scaling results above, we conclude that the nonlinear observations here result from a direct intensity or momentum fluctuation, rather than a secondary instability arising from a dark stripe.
Compressibility plays an interesting role in these experiments. Its influence can be seen by comparing the growth rate of a sinusoidal interface perturbation with wavelength l, ω g ∼ (g/l) 1/2 , to the relaxation rate of density waves at the same wavelength ω s ∼ c s /l, where c s ∼ 2 µm mm −1 denotes the effective sound speed. Simple scaling suggests that compressibility effects should become dominant when ω s ∼ ω g , or l ∼ l * c 2 s /g. For the system here, l * is of the order of 4 µm, while the measured perturbation periods are of the order of 100 µm, implying that the instability dynamics should be fully compressible. As shown by equation (5) and figure 4, however, compressibility appears dominant only in one direction, normal to the interface. Physically, this follows by realizing that shock dynamics requires an initial density difference [16]; for the basic RTI geometry, such a jump occurs across the interface but not along it.
We emphasize that compressibility is a singular perturbation, involving directly the equation of state (and thus nonlinearity), and that even in normal fluids, including compressibility is a difficult and somewhat controversial task [48]. The consensus is that compressibility suppresses the RTI, as energy that would drive perturbations is now redirected into squeezing the fluid [49,50]. However, other arguments suggest that the finite sound speed from compressibility causes energy to accumulate near the interface (rather than radiate from it), thus destabilizing the flow [51,52]. Here, the observations suggest that compressibility is a destabilizing influence in quantum RTI, as the initial shock wave enables the interface to travel faster than it would without the shock. As this transport occurs in the same direction as 'gravity', the end result is a higher effective acceleration, resulting in a higher-energy perturbation. While the gain formula (5) captures the essence of these results, simply by substituting the appropriate shock length, a more quantitative explanation of our experimental observations can only be provided by a more complete theory of compressible Schrödinger flow, which is yet to be developed.
Finally, we note that the general form of the equations and experimental setup means that the approach here can be extended in many possible ways, including flow with multiple layers, fully three-dimensional dynamics and partially condensed (statistical) systems. It also suggests that light-matter interactions in more general environments, such as optofluidic systems, can be treated using a unified language in which the fields are regarded as equal components in a multispecies fluid, rather than as separate entities.

Conclusions
We have used a coherent optical system to experimentally explore the dynamics of the superfluid RTI. Parametric studies were conducted with density difference, acceleration and nonlinear interaction strength, revealing behavior that was compressible across the interface but effectively incompressible parallel to it. In particular, data collapse occurred when shock-like scaling for the mixing layer was substituted into the incompressible gain formula. The results hold for any Schrödinger fluid, e.g. superfluids and quantum plasma, and lay the foundation for a variety of fluid-inspired instabilities in optics.