Turbulence in a matter-wave supersolid

Quantum turbulence associated with wave and vortex dynamics is numerically investigated for a two-dimensional trapped atomic Rydberg-dressed Bose-Einstein condensate (BEC). When the coupling constant of the soft-core interaction is over a critical value, the superfluid (SF) system can transition into a hexagonal supersolid (SS) state. Based on the Gross-Pitaevskii equation approach, we have discovered a new characteristic k−13/3 scaling law for wave turbulence in the SS state, that coexists with the waveaction k−1/3 and energy k−1 cascades commonly existing in a SF BEC. The new k−13/3 scaling law implies that the SS system exhibits a negative, minus-one power energy dispersion (E ~ k−1) at the wavevector consistent with the radius of the SS droplet. For vortex turbulence, in addition to the presence of the Kolmogorov energy k−5/3 and Saffman enstrophy k−4 cascades, it is found that large amount of independent vortices and antivortices pinned to the interior of the oscillating SS results in a strong k−1 scaling at the wavevector consistent with the SS lattice constant.

In addition to VT, turbulence consisting of waves, named wave turbulence (WT), also plays a central role in a system of interactions. As a general phenomenon, WT is observed in a vast range of nonlinear systems, on quantum to astrophysical scales 28 . In a 2D atomic SF BEC 29 , WT scaling laws typically involve the k −1/3 waveaction-cascade and the k −1 compressible energy-cascade spectra 29 . The former corresponds to the strongly nonequilibrated process of evaporative cooling and the latter corresponds to the condensation process, respectively. Supersolid (SS) is a state of matter that simultaneously possesses superfluidity and solidity and in which both gauge and continuous translational symmetries are broken [30][31][32][33][34] . Recently SS states have been realized in a BEC coupled to the modes of two optical cavities 35 and in a BEC with spin-orbit coupling 36 . Another promising candidate for the SS state is the Rydberg-dressed atomic BEC that exhibits a defocusing soft-core interaction [37][38][39][40][41] . In such a system, the interaction between the Rydberg-dressed ground-state atoms behaves as: the defocusing interaction strength, R c the blockade radius, and ′ ≡ − r r r the relative position of two particles. Let Λ =  C 6  denote the coupling constant and when Λ is over a critical value Λ c , the 2D SF system can transition into a hexagonal SS state. Figure 1 schematically shows the ground states of the Rydberg-dressed BEC with the coupling constant below (SF state) and above (SS state) the critical value. Compared to the SF state, there are two new emerging length scales in the SS state: the lattice constant, d, and the radius of the SS droplets, R, which play important roles in QT of a SS.
The main theme of this paper is to investigate how the WT and VT behave in a SS state. Are there any new scaling laws characterizing the SS state? What roles are the new length scales playing in such a state? It is worth noting that, a self-organized structure associated with density modulation also exist in a variety of systems from Here V pot (r) = mω 2 r 2 /2 is the harmonic trapping potential with ω the frequency and m the atom mass and as introduced earlier, 6 is the defocusing soft-core interaction. The order parameter ψ, satisfying the normalization condition ∫ ψ = dr 1 2 , is the condensate wave function. Throughout this paper, the blockade radius R c and are used as the units of length and time, respectively. In our simulation, a SF condensate with Thomas-Fermi (TF) radius R TF = 6R c was initially prepared. For the SS ground state, coupling constant is chosen to be Λ = 2400, above the critical value Λ c = 2000. The trapping frequency is set at ω = 3. To generate turbulence, circulation-3 vortex and antivortex were imprinted. (Note that more initial vortex pairs have been tested and essentially lead to the same results.) Time evolution of the wave function follows the GPE: iℏ∂ψ/∂t = δE/δψ * . To numerically integrate GPE, we use the method of lines with spatial discretization by highly accurate Fourier pseudospectral method and the time integration is done by adaptive Runge-Kutta method of orders 4 and 5 (RK45). As resolving the trap potential can reduce the available numerical power, our simulations are performed on a square grid of 512 2 points in a 30 2 domain.
About the error of simulating GPE, it can be discussed by spatial and time errors of our numerical method. For spatial accuracy, we have employed Fourier-pseudospectral method, which has the property of exponential convergence compared with traditional algebraic-order-convergence methods like finite difference or finite element. High-order pseudospectral methods generally provide excellent spatial accuracy with economically practicable resolutions. For time integration, we used RK45, a time-step-adjustable Runge-Kutta method to meet the specified error tolerance, which automatically satisfies the stability criterion as well. Under these treatments and some benchmark validation, we are assured that the wave function has not been contaminated by numerical noise (round-off error) even after long-time computation.
We shall focus on the turbulence shown in the kinetic energy spectrum, E kin (t). As to extract the results for both WT and VT, one can separate E kin (t) into two mutually orthogonal parts, kin , i.e., compressible and incompressible. In this regard, it is useful to express the condensate wave function ψ(r, t) in can then be decomposed into irrotational and solenoidal parts, or correspondingly, compressible and incompressible parts [47][48][49][50] Because n u ( ) c and n u ( ) i are mutually orthogonal, the decomposition of the kinetic energy density follows: , 2 . Physically  c kin and  i kin correspond respectively to the kinetic energy densities of the sound wave and the swirls in a superflow.
To study the scaling laws of the turbulence, it is to to the k space using the sum rule is defined as the angle-averaged kinetic-energy spectrum  Wave turbulence. In view of Fig. 2(a) for the changes in the compressible energy, we have identified that a waveaction-cascade k −1/3 spectrum occurs at k s < k < k M , which corresponds to the condensation process. An energy-cascade k −1 spectrum at k M < k < k R is also identified, which corresponds to energy transport away from the interior of the condensate through splitting into small-scale sound waves. The lower limit of the k −1/3 waveaction-cascade spectrum is k s = 2π/s ≈ 1.8, where s ≈ 3.5 corresponds to the eventual condensate size and represents the longest length scale of the system. The upper limit of the k −1/3 scaling law is near the wavevector k M . As shown in Fig. 1 in the firat Brillouin zone of a 2D hexagonal lattice,

Scaling Laws
. The lower (upper) limit of the k −1 scaling law is near the wavevector k M (K R ), where k R = 2π/R ≈ 21.5 with R corresponding to the radius of the SS droplet (see Fig. 1). The k −1/3 waveaction-cascade and k −1 energy-cascade are the two scaling laws commonly seen in a 2D SF BEC. Both spectra correspond to a quadratic energy dispersion E ~ k 2 at low k in an N = 4 process 28 .
In addition to the above two scaling laws, a new k −13/3 scaling law is seen to appear uniquely for the SS state in the ultraviolet regions, starting at k R . By examining and comparing the compressible energy spectrum near t = 75 (before equilibrium) to that near t = 150 (after equilibrium), it is found that the ultraviolet energy spectrum increases at smaller k's and decreases at larger k's in the time frame. It means that the direction of energy flux in k space is from larger k to smaller k 28 . We thus conclude that this new −13/3 scaling corresponds to an inverse waveaction-cascade. We propose that the new −13/3 scaling law is the signature in WT for a 2D SS. The lower limit of this spectrum at k R is consistent with the radius of the SS droplet. The upper limit is expected to be bounded by the shortest length scale of the system, i.e., the healing length ξ ≈ 0.05. It will correspond to a wavevector k ξ = 2π/ξ ≈ 125. In Fig. 2(a), the calculated −13/3 spectrum for a 2D trapped Rydberg-dressed BEC is only shown for wavevectors up to k ~ 35. Due to the effect of computing resolution, the spectrum deviates from a straight line at  k 50. To verify whether this new scaling law can sustain a larger k range when computing resolution is increased, we perform a calculation on a similar system without the trap (by setting trapping frequency ω → 0). The kind of calculation enables one to simulate an infinite system by using a unit cell. In this case, the result of WT spectrum is shown in Fig. 3 for the ultraviolet regime. The k −13/3 spectrum is seen to cover almost a decade, ranging from  k 15 to 80. The compressible and incompressible sectors can possibly interact in the presence of the static spatial modulation. Therefore it is important to justify the new −13/3 scaling law as the value is close to the Saffman −4 law appearing in the incompressible sector [see Fig. 2(b)]. We have applied the method of least squares within the relevant scales and the line of best fits, to the results in Fig. 3, appears to have a slope −4.318. This number is consistent with the proposed −13/3. In any case, one should be more careful. The value −4.318 quoted implies the accuracy in the fourth digit and the exponent −13/3 coincides with −4.318 only up to the second digit, thus the error bar on the exponent is about 0.01 (or 2%). However, if taking into account the short-scale fluctuations happening near k = 60 (see Fig. 3) which covers roughly 1/4 of the order, and when comparing to best-fit line that covers about 4 orders of magnitude, the relative error of the exponent could be up to 1/16 (or 6%). Thus the error bar on the −13/3 law could be up to ~0. 25. This means that the law −13/3 can well be −4.0 or −4.5 and accordingly, the appearance of the Saffman −4 law can not be excluded.
The new k −13/3 waveaction-cascade associated with the formation of SS is very distinct from the k −1/3 waveaction-cascade associated with the formation of condensate. The lower bound of the new k −13/3 scaling coincides with the radius R of the SS droplet, while the lower bound of the k −1/3 scaling coincides with the eventual radius s of the condensate. Vortex turbulence. In the incompressible VT energy spectra displayed in Fig. 2(b), both Kolmogorov k −5/3 and Saffman k −4 scaling laws are identified. The Kolmogorov k −5/3 scaling law occurs at k d′ < k < k a . The lower limit k d′ = 2π/d′ ≈ 8 with d′ = d − 2R ≈ 0.8 corresponding to the gap between two adjacent SS droplets (note that when R → 0, d′ → d). Whereas the upper limit k a = 2π/a ≈ 14 with a ≈ 0.45 corresponding to the radius of the maximum circle formed in the interior of the three adjacent SS droplets. The geometry gives a simple relationship between the three characteristic lengths: The radius a also corresponds to the range over which the fluctuating vortices are mainly located [see Fig. 4(c,d)]. The lower limit of the Saffman k −4 scaling law occurs at k a , while the upper limit is expected to be bound by the largest k scale, k ξ , as well.
In addition to the Kolmogorov k −5/3 and Saffman k −4 scaling laws commonly seen in VT for a 2D SF BEC, we also identify a strong k −1 scaling law that covers almost a decade of the k range in the infrared regime. The k −1 spectrum corresponds to the motions of many independent vortices or antivortices 27,51,52 . Oscillation of the SS can constantly create vortices and antivortices in the interior, and such a fast vortex-antivortex creation and annihilation cycle results in the k −1 scaling. The lower limit of the k −1 scaling also coincides with the longest length scale s of the system, while the upper limit coincides with k d′ .
Because hexagonal lattice is one of the most critical characteristics of the present system, lattice constant d or correspondingly wavevector k d plays a crucial role in both the compressible and incompressible spectra. One sees in Fig. 2(a,b) that energy peaks rise and fall at k = k d in both the compressible and incompressible spectra. The rise and fall of peaks are especially obvious in the incompressible one, which correspond to the creation and annihilation of vortices in the interior of the SS structure.
Local Energy density distribution. It is also important and interesting to see how the order parameter and energy density distribute in real space. At time t = 150 after the equilibrium, Fig. 4(a) plots the normalized wave function (the order parameter) |Ψ(r)| in a square space. Both the initial and eventual sizes (R TF = 6R c and s ≈ 3.5R c ) of the condensate are displayed. SS structure is seen to be robustly sustained. Figure 4(b,c) present the corresponding compressible and incompressible local energy densities,  r ( ) c kin and  r ( ) i kin . The compressible and incompressible local energy distributions are seen to overlap strongly not only in the border area but also in the interior. The inset in Fig. 4(b) shows clearly that the radius R of the SS droplet plays an important role in the WT of a SS. This also explains why the lower bound of the new k −13/3 scaling law coincides with k R .
Because of the formation of the hexagonal SS lattice, the fluctuating vortices shown in Fig. 4(c) for VT are seen to form an analogous hexagonal structure in the interior part. That is, the density distribution of the vortex core is highly inhomogeneous. This explains why the lower bound of the Saffman k −4 scaling coincides with a short length scale a. To see this in more details, we plot in Fig. 4(d) the stream function φ(r) associated with the velocity field n u 53 . It is the solution of Poisson's equation φ ∇ = ⋅ Ω → z 2 , where Ω → = ∇ × n u ( ) is the 2D vorticity vector. In Fig. 4(d), the colors indicate the corresponding values of φ(r), and the locations and structures of eddies can be easily recognized by families of closed streamlines. From the definition of the vorticity vector Ω → , we note that the singular and regular parts of Ω → are modulated by n and its gradient, respectively. As a result, the magnitude of the vorticity Ω → of a vortex in high-density area is greater than that of a vortex in low-density area and, with such reasoning, we state that an eddy is energetically "larger" if it contains more circumfluent particles. Figure 4(d) clearly shows that the streams are mostly concentrated in a circle of radius a, which indicates that the fluctuating vortices mainly reside in this region. The figure also reveals that the streamline bundles form a (fluctuating) hexagonal structure, which is consistent with the incompressible VT energy spectra displayed in Fig. 4(c).

On the −13/3 scaling law
It is important to investigate what causes the new WT −13/3 scaling law? As mentioned before, the −13/3 scaling law corresponds to an inverse waveaction cascade, so we proceed to search for what wave is responsible for. Let us first consider a wave that has an effective dispersion ω k ~ λk α for certain k range, where λ is a positive parameter and α is the power. The dimension of λ is thus  28 . In an N-wave process, one also has the rate relation, ε ∼ ∼ −  E E N 1 . In addition to the conserved energy, for an even number N-wave process (as for the GPE simulation, N = 4), the total waveaction or "particle number", ∫ ∫ , is also conserved for the N/2 → N/2 processes. Thus one can also have the following waveaction balance equation, ∂ t n k + ∂ k ζ = 0, where ζ = ε/ω k is the waveac- It is of great interest if one can solve the entire wave dispersion for the current Rydberg-dressed BEC system and confirm that ω k ~ λk −1 for k ~ k R . This is not feasible, however. We instead seek for ideas from the elementary excitation to which analytical results are available for the system 32 . In the absence of the trapping potential, the stable elementary excitation of the system in the SF state (when coupling constant Λ < Λ c ) behaves as where n 0 is the uniform density and = | | ∼ ∼ U k U k ( ) ( ) is the Fourier transformation of the isotropic soft-core interaction function U(r) introduced earlier. Since ∼ U k ( ) is negative in certain k intervals, dispersion (5) could become unstable when Λ > Λ c . It means that above the critical point, an instability, named "roton instability", occurs and the SF state will transition into a periodic SS state. In Fig. 5, based on Eq. (5) the blue and red curves plot the elementary excitation spectra ω(k) for the stable SF state at Λ Λ  c and  Λ Λ c , respectively. Whereas orange curve plots the longitudinal elementary excitation spectrum ω(k) for the periodic hexagonal SS state at Λ > Λ c which is obtained by solving the Bogoliubov-de Gennes (BdG) equation. Of most interest, a low-energy gapless mode is seen to appear at → − k k M for the SS state, where k M is defined in Fig. 1. As displayed by the black line in Fig. 5, the asymptotic behavior of the orange curve at → − k k M does exhibit a minus-one power dispersion, ω(k) ~ k −1 .
The low-energy gapless k −1 elementary excitation does not have a direct relation to the new −13/3 WT scaling law though. One reason is that the low-energy gapless k −1 elementary excitation occurs at k M , while the k −1 wave dispersion that causes the new −13/3 WT scaling law occurs at k R . Secondly, the wave that causes the new −13/3 WT scaling law should have relatively high energy, while the k −1 elementary excitation is the low-energy one. Nevertheless, it still shed light on why in the SS state, the wave dispersion behaves ~k −1 in the ultraviolet regions. When Λ > Λ c in the SS state, negative part of ∼ U k ( ) can result in a significant reduction in the energy spectrum. The free-particle wave dispersion k 2 for the kinetic energy can be significantly reduced and consequently a down-turn ω(k) can occur. A down-turn dispersion implies a negative power. As for why the new −13/3 WT scaling law coincides with k R , it can be understood in this way. When SS forms, the condensed atoms pile up in a lattice of droplets. Thus the ultraviolet waves will conform to the droplets with an onset wavevector k R .

Conclusion
In summary, quantum turbulence including both wave and vortex dynamics is investigated in a 2D trapped Rydberg-dressed Supersolid (SS). The SS state of the Rydberg-dressed BEC occurs when the coupling constant is over a critical point. Of most interest, we discover a new −13/3 scaling law in the wave turbulence, which is the signature for the 2D SS state. This −13/3 scaling implies that the wave has an effective k −1 dispersion at the wavevector consistent with the radius of the SS droplets. As self-organized structure associated with density modulation can exist in a variety of systems from nonlinear optics 42 , water waves 43 , plasma 44,45 , to planetary waves 46 , our study should shed some light on the wave turbulence of other systems.