Kolmogorov-Hinze scales in turbulent superfluids
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] J a n Kolmogorov-Hinze scales in turbulent superfluids
Tsuyoshi Kadokura and Hiroki Saito Department of Engineering Science, University of Electro-Communications, Tokyo 182-8585, Japan (Dated: January 7, 2021)When a two-component mixture of immiscible fluids is stirred, the fluids are split into smallerdomains with more vigorous stirring. We numerically investigate the sizes of such domains in astirred two-component superfluid in a fully-developed turbulent state with energy-input rate ǫ . Forthe strongly immiscible condition, the typical domain size is shown to be proportional to ǫ − / , aspredicted by the Kolmogorov-Hinze theory in classical fluids. For the weakly immiscible condition,quantum effects become pronounced and the power changes from − / − /
4. The domain sizescaled by the interface thickness collapses into a universal line.
When oil is poured into water and these fluids arestirred, the oil becomes split into droplets in the water.The droplet sizes become smaller with more vigorous stir-ring. Such disintegration phenomena in multicomponentfluids are ubiquitous in nature, and are important in en-gineering and industry.Kolmogorov [1] and Hinze [2] considered the disinte-gration process of droplets, and estimated the size ofdroplets in turbulent fluids. In fully-developed turbu-lence, the energy is input into the system as large-scaleeddies that cascade toward a smaller scale, resulting inthe Kolmogorov power law of the energy spectrum [3].In such turbulent fluids, large-size droplets are unstablebecause they are fragile to deformation and disintegra-tion due to the fluctuating pressure of the surroundingfluid. Small droplets are thus produced by the breakupof large droplets, and this breakup process continues to ascale where the turbulent energy to break up the dropletsbecomes balanced with the droplet energy sustaining itsshape. Droplets smaller than this scale coalesce into largedroplets. Therefore, there exists a characteristic size D of droplets in turbulent fluids, which is referred to as theKolmogorov-Hinze (KH) scale, given by [2] D ∼ ( σ/ρ ) / ǫ − / , (1)where σ is the interface tension coefficient, ρ is the den-sity of the surrounding fluid, and ǫ is the energy inputrate to maintain the turbulence. The KH scale has beenexperimentally verified in various systems [4–8]. Further-more, direct numerical simulations have been performedover the last decade [9–14].In this Letter, we extend the study of KH scales to aquantum mechanical system: the superfluid turbulence ofa two-component Bose-Einstein condensate (BEC). Wewill show that the KH scale also appears in this superfluidsystem and is modified by quantum effects. Turbulent be-havior in superfluids has been widely studied. For single-component superfluids, an isotropic fully-developed tur-bulent state exhibits the Kolmogorov power law [15–20].The turbulent behavior of gaseous BECs has also beenexperimentally studied [21–23], and a power law behav-ior has been observed recently [24–27]. A wide varietyof systems have been studied theoretically, such as two-dimensional systems [28–32], dipolar superfluids [33], and boundary layers [34]. Here, we focus on the turbulencein a two-component BEC. Turbulence in multicomponentBECs has been investigated by many researchers [35–41].In the context of domain-size scaling in multicomponentBECs, coarsening dynamics following domain formationhave been studied extensively [42–52]. However, the KHscale, i.e., domain-size scaling in conjunction with theKolmogorov turbulence, has not yet been investigated.The KH scale in Eq. (1) is derived as follows. In a tur-bulent fluid, a domain undergoes fluctuating pressuresthat vary over its size D , which causes deformation anddisintegration of the domain. This pressure differenceover the size D can be expressed as ∼ ρv ≡ P turbulence ,where v is the velocity difference of the surrounding fluidover the size D . Within the inertial range of an isotropichomogeneous turbulence, the statistical average of v obeys the Kolmogorov two-thirds law [3], ¯ v ∝ ( ǫD ) / ,and hence P turbulence ∼ ρ ( ǫD ) / . On the other hand,a domain tends to sustain its shape and resist disinte-gration. This sustaining force arises from the interfacetension, and the pressure required to deform the domainis thus estimated to be ∼ σ/D ≡ P sustain [53]. Breakupof domains to smaller sizes stops at the scale that satisfies P sustain ∼ P turbulence , (2)which gives the KH scale in Eq. (1).For an immiscible two-component BEC, the interfacetension, which arises from the interatomic interactionand quantum pressure, is well-defined, as in classical flu-ids [54–56]. Therefore, we expect that the KH scale inEq. (1) also emerges in two-component BECs when thethickness of the interface W is much smaller than thedomain size D . However, when W is comparable to orlarger than D , the picture of the interface tension breaksdown in the derivation of Eq. (1). The interface thick-ness W is determined by the competition between thequantum pressure and the intercomponent repulsion, and W becomes large when the former dominates the lat-ter. Thus, in the limit of large W , the mechanism tosustain domains against disintegration originates mainlyfrom the quantum pressure P sustain ∼ ~ / ( mD ) insteadof P sustain ∼ σ/D , where m is the mass of an atom. Inthis case, Eq. (2) gives the characteristic size as D ∼ ( ~ /m ) / ǫ − / . (3)Therefore, in the limit of weak segregation with large W ,the quantum mechanical effect becomes pronounced andthe KH scale is expected to change from the − / − / ǫ . In the remainder ofthis Letter, we will corroborate this prediction using nu-merical simulations of the coupled Gross-Pitaevskii (GP)equations.In the mean-field approximation, a two-componentBEC is described by the coupled GP equations,( i − γ ) ~ ∂ψ ∂t = (cid:18) − ~ m ∇ + V ext + g | ψ | + g | ψ | − µ (cid:19) ψ , (4a)( i − γ ) ~ ∂ψ ∂t = (cid:18) − ~ m ∇ + V ext + g | ψ | + g | ψ | − µ (cid:19) ψ , (4b)where ψ j ( r , t ) is the macroscopic wave function of the j th component, V ext ( r , t ) is the external stirring poten-tial, and g jj ′ = 4 π ~ a jj ′ /m with a jj ′ being the s -wavescattering length between the j th and j ′ th components.The constant γ in Eq. (4) introduces dissipation of energyfrom the condensate into thermal atoms [57]. The valueof γ is selected in such a way that the energy dissipa-tion occurs predominantly below the scale of the healinglength and does not affect the dynamics in the inertialrange. The chemical potential µ j is adjusted to maintainthe unitarity, i.e., R | ψ j | d r is conserved.The miscibility between the two components is deter-mined by the coupling coefficients g jj ′ . The two super-fluids are immiscible and phase separation occurs when g > g g is satisfied [58]. In the following, we assume g = g ≡ g > g >
0; therefore, the immisciblecondition reduces to g > g. (5)The phase separation of immiscible components producesan interface, at which excess energy arises, resulting inthe interface tension. For g /g − ≪
1, the interfacetension coefficient is given by [54–56] σ ≃ (cid:20) ~ n m ( g − g ) (cid:21) / , (6)where n is the number density far from the interface. Thethickness W of the interface between two components,over which each density changes from 0 to n (or n to 0),has the form W ≃ ξ ( g /g − − / , (7)where ξ = ~ / (2 mgn ) / is the healing length.The length, time, and wave functions are normalizedby ξ , ~ / ( gn ), and √ n , respectively, where n is takento be an average density of each component. In this FIG. 1. A two-component BEC is stirred by a plate-shapedpotential (red or black) to generate a turbulent state. (a)Snapshot of a fully-developed turbulent state at t = 2000 for g ′ ≡ g /g = 1 .
02 and ν plate = 0 . | ψ | = 0 . | ψ | = 0 . | ψ | and | ψ | + | ψ | at z = 32. The arrows indi-cate the direction of the potential motion. (d) Energy spectra E ( k ) (arbitrary unit) for γ = 0 . , .
03, and 0 .
05, averagedover t = 1900-2000. The dashed line has a slope of − /
3. Seethe Supplemental Material for a movie of the dynamics [60]. unit, the interaction coefficients appear in the normal-ized GP equation only through g /g ≡ g ′ . We considera box of size L = 64 with a periodic boundary condi-tion, in which the two components are equally populated, R | ψ | d r = R | ψ | d r = L . To generate the turbulentstate, the system is stirred by a plate-shaped potential V ext with a size of 1 × ×
64 and a potential heightof 10. This plate-shaped potential is moved sinusoidallyover the system (0 < x < L ) at a frequency ν plate , asshown in Fig. 1(a), which can input the turbulent energyinto the large scale. The coupled GP equations (4) arenumerically solved using the pseudospectral method [59].The space is discretized into a 128 mesh and the timestep is typically 10 − . The initial state has uniform den-sity with random phases on each mesh.The two-component BEC is stirred until the fully-developed turbulent state is achieved. Figure 1(a) showsisodensity surfaces of | ψ | and | ψ | at t = 2000 for g ′ = 1 .
02 and ν plate = 0 . g ′ ; therefore, thetwo components are separated and domains are formed ineach component. Figure 1(a) shows that the domain sizesare typically > ∼
10, while W ≃ (1 . − − / ≃ | ψ | and | ψ | + | ψ | , respectively, where the stirring poten-tial is moving rightward at this moment. Although thedensity | ψ | (or | ψ | ) largely varies in space due to thephase separation (Fig. 1(b)), the total density far fromthe stirring potential is almost uniform (Fig. 1(c)). Thesmall density holes in Fig. 1(c) indicate that quantizedvortices are generated behind the potential moving atalmost the sound velocity of the density wave. To con-firm that the system has reached the Kolmogorov turbu-lence, we calculate the energy spectrum E ( k ) of incom-pressible velocity field [15] averaged over the duration of t = 1900-2000 for γ = 0 .
01, 0.03, and 0.05, which isshown in Fig. 1(d). The slope of the logarithmic plot of E ( k ) agrees well with the Kolmogorov power law of − / k < ∼ γ , which indicates that, forthese values of γ , the energy dissipation occurs predom-inantly at the small length scale and does not affect thedynamics in the inertial range. In the following calcula-tions, we take γ = 0 . d ( r ) = (cid:26) | ψ ( r ) | ≥ | ψ ( r ) | ) − | ψ ( r ) | < | ψ ( r ) | ) , (8)which eliminates the information of the interface thick-ness and small excitations, and enables the informationof only the domain size to be extracted. Figures 2(a)-2(c) show snapshots of the distributions of d ( r ). Fouriertransformation of d ( r ) is performed and its azimuthalaverage in the k space is taken, which we denote as˜ d ( k ). Figure 2(d) shows the time evolution of ˜ d ( k ) forthe dynamics in Fig. 1. Since the initial state has a ran-dom phase on each mesh, ˜ d ( k ) is initially independentof k . The distribution of ˜ d ( k ) then tilts toward small FIG. 2. Discretized density distributions d ( r ) in Eq. (8)at (a) t = 1 for g ′ ≡ g /g = 1 .
02, (b) t = 2000 for g ′ = 1 .
02, and (c) t = 2000 for g ′ = 1 . d ( k ) for g ′ = 1 .
02. In (a)-(d), ν plate = 0 . D ≡ π/ [ R k ˜ d ( k ) dk/ R ˜ d ( k ) dk ] for( g ′ , ν plate ) = (1 . , . . , . . , . . , . wave numbers, and the steady distribution is reached for t > ∼ d ( k ) does not depend on the initial state. Thefirst moment of ˜ d ( k ), R k ˜ d ( k ) dk/ R ˜ d ( k ) dk ≡ h k i , givesa typical wave number in the density distribution, andthus 2 π/ h k i ≡ D characterizes the size of the domains.Figure 2(e) shows the time development of D for differ-ent values of g ′ and ν plate . The size D in the steady statedecreases with g ′ and with an increase in ν plate . This isbecause the interface tension (6) decreases with g ′ andthe energy input rate increases with ν plate .To study the KH scale quantitatively, we need to cal-culate the rate ǫ at which the energy is input into thesystem by the stirring potential. Let the position ofthe plate-shaped potential be x plate ( t ), which moves inthe ± x direction, as shown in Fig. 1. The force act-ing on the potential from the superfluid is then F = − ∂/ ( ∂x plate ) R V ext ( r ; x plate )[ | ψ ( r ) | + | ψ ( r ) | ] d r . Theenergy input rate is thus given by ǫ = − F ˙ x plate . Thevalue of ǫ (and also D ) fluctuates in time due to the ran-dom nature of the turbulence. The sinusoidal motionof the plate-shaped potential also causes periodic fluc-tuation, as shown in Fig. 2(e). Therefore, we take thetemporal average of these quantities, ¯ ǫ and ¯ D , over asufficiently long time after the steady turbulent state isachieved. The inset in Fig. 3(a) shows ¯ ǫ as a function of ν plate for g ′ = 1 .
02 and 1.001.
FIG. 3. (a) Typical domain size ¯ D versus energy-inputrate ¯ ǫ for the steady turbulent state with g ′ ≡ g /g = 1 . − / − / ǫ as afunction of ν plate . (b) ¯ D and ¯ ǫ are rescaled with respect tothe interface width W for various values of g ′ . The plotscollapse into a universal curve. The slopes of the lines are − / − / Now we are ready to investigate the KH scales in a tur-bulent superfluid in the classical and quantum regimes,as given in Eqs. (1) and (3). The result is shown in Fig. 3,which is the main result of this work. Figure 3(a) plotsthe typical domain size ¯ D versus the energy-input rate ¯ ǫ for g ′ = 1 .
02 and 1.001. For g ′ = 1 .
02, the plots obey thepower law ∝ ¯ ǫ − / , which agrees with the classical KHscale in Eq. (1). This implies that the two componentsare well separated, and the mechanism that sustains thedomains against disintegration can be described by theinterface tension in this region. For g ′ = 1 . ∝ ¯ ǫ − / , which agrees with the KH scale in the quantumregion in Eq. (3), implying that the mechanism to sustaindomains is mainly the quantum kinetic pressure arisingfrom the uncertainty principle. In this region, ¯ D < ∼
10 issmaller than W ≃ (1 . − − / ≃
32 and the interfacetension of domains is not well-defined.As discussed in deriving Eqs. (1) and (3), and numer- ically confirmed in Fig. 3(a), the different power lawsemerge, depending on the two length scales: the domainsize ¯ D and the interface width W = ( g ′ − − / . There-fore, we use the rescaled domain size, ¯ D/W . We candeduce that the crossover between the two power laws,Eqs. (1) and (3), always lies in the region of ¯
D/W ∼ W , and therefore,in the quantum KH region, ¯ D = c ¯ ǫ − / can be rewrit-ten as ¯ D/W = c ( W ¯ ǫ ) − / with a coefficient c that isindependent of W . The ¯ D/W versus W ¯ ǫ line is thusuniversal in the region of the − / − / D/W ∼
1; therefore it is expected that this line isuniversal over the entire region. Figure 3(b) plots ¯
D/W versus W ¯ ǫ for various values of g ′ . As expected, all theplots collapse into the universal curve, which bridges thetwo regions of the − / − / D and ¯ ǫ are restricted to narrow ranges for each valueof g ′ in the present numerical simulations, the rescaledplot in Fig. 3(b) significantly extends the effective rangesof ¯ D and ¯ ǫ , which corroborates the existence of the twopower laws in the superfluid KH scale.Finally, we discuss the possible experimental realiza-tion of the present result. A box potential will be suit-able to avoid complexity arising from the inhomogeneous | ψ | + | ψ | distribution in a harmonic potential. Thestirring potential can be produced by a far-off-resonantlaser beam. Shaking of an optical box can also be used togenerate the turbulent state [25]. The typical size of thedomains can be inferred from the imaging data, where theslice imaging of a three-dimensional distribution may berequired [61]. It is difficult to measure the energy inputrate directly; therefore, the support of numerical simula-tion is necessary, which provides the relation between themotion of the potential and the energy input rate, as inthe inset in Fig. 3(a). The interaction g can be variedusing the Feshbach resonance technique.In conclusion, we have investigated the KH scale ofdomain sizes in immiscible two-component superfluids ina fully-developed turbulent state. We predict that tworegions of the KH scale exist with different power laws,which reflect the quantum property of the system. Nu-merical simulations of the coupled GP equations wereperformed, and the typical domain size D was confirmedto obey the power laws with respect to the energy inputrate ǫ . We found the universal function of the rescaledquantities D/W and W ǫ , which connects the two re-gions of the quantum and classical KH scales. A possibleextension of this study is a three-component system, inwhich the third component can change the interface ten-sion of the other two components, resulting in emulsifi-cation.This work was supported by JSPS KAKENHI GrantNumber JP20K03804. [1] A. N. Kolmogorov, Dokl. Akad. Nauk. SSSR , 825(1949). [2] J. O. Hinze, AIChE J. , 289 (1955). [3] See, e.g., U. Frisch, Turbulence: The Legacy of A. N. Kol-mogorov (Cambridge University Press, Cambridge, UK,1995).[4] P. H. Clay, Proc. Roy. Acad. Sci. (Amsterdam) , 852(1940).[5] R. Shinnar, J. Fluid Mech. , 259 (1961).[6] C. A. Sleicher, AIChE J. , 471 (1964).[7] K. Arai, M. Konno, Y. Matunaga, and S. Saito, J. Chem.Eng. Jpn. , 325 (1977).[8] G. B. Deane and M. D. Stokes, Nature (London) ,839 (2002).[9] P. Perlekar, L. Biferale, M. Sbragaglia, S. Srivastava, andF. Toschi, Phys. Fluids , 065101 (2012).[10] R. Skartlien, E. Sollum, and H. Schumann, J. Chem.Phys. , 174901 (2013).[11] P. Perlekar, R. Benzi, H. J. H. Clercx, D. R. Nelson, andF. Toschi, Phys. Rev. Lett. , 014502 (2014).[12] X. Fan, P. H. Diamond, L. Chac´on, and H. Li, Phys. Rev.Fluids , 054403 (2016).[13] P. Perlekar, N. Pal, and R. Pandit, Sci. Rep. , 44589(2017).[14] M. E. Rosti, Z. Ge, S. S. Jain, M. S. Dodd, and L. Brandt,J. Fluid Mech. , 962 (2019).[15] C. Nore, M. Abid, and M. E. Brachet, Phys. Rev. Lett. , 3896 (1997); Phys. Fluids , 2644 (1997).[16] S. R. Stalp, L. Skrbek, and R. J. Donnelly, Phys. Rev.Lett. , 4831 (1999).[17] T. Araki, M. Tsubota, and S. K. Nemirovskii, Phys. Rev.Lett. , 145301 (2002).[18] M. Kobayashi and M. Tsubota, Phys. Rev. Lett. ,065302 (2005).[19] N. G. Parker and C. S. Adams, Phys. Rev. Lett. ,145301 (2005).[20] A. W. Baggaley, J. Laurie, and C. F. Barenghi, Phys.Rev. Lett. , 205304 (2012).[21] E. A. L. Henn, J. A. Seman, G. Roati, K. M. F. Mag-alh˜aes, and V. S. Bagnato, Phys. Rev. Lett. , 045301(2009).[22] T. W. Neely, A. S. Bradley, E. C. Samson, S. J. Rooney,E. M. Wright, K. J. H. Law, R. Carretero-Gonz´alez, P.G. Kevrekidis, M. J. Davis, and B. P. Anderson, Phys.Rev. Lett. , 235301 (2013).[23] W. J. Kwon, G. Moon, J.-Y. Choi, S. W. Seo, Y.-I. Shin,Phys. Rev. A , 063627 (2014).[24] K. J. Thompson, G. G. Bagnato, G. D. Telles, M. A.Caracanhas, F. E. A. dos Santos, and V. S. Bagnato,Laser Phys. Lett. , 015501 (2014).[25] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic,Nature (London) , 72 (2016).[26] S. P. Johnstone, A. J. Groszek, P. T. Starkey, C. J.Billington, T. P. Simula, and K. Helmerson, Science ,1267 (2019).[27] N. Navon, C. Eigen, J. Zhang, R. Lopes, A. L. Gaunt, K.Fujimoto, M. Tsubota, R. P. Smith, and Z. Hadzibabic,Science , 382 (2019).[28] S. Nazarenko and M. Onorato, J. Low Temp. Phys. ,31 (2007).[29] T.-L. Horng, C.-H. Hsueh, S.-W. Su, Y.-M. Kao, andS.-C. Gou, Phys. Rev. A , 023618 (2009).[30] R. Numasato, M. Tsubota, and V. S. L’vov, Phys. Rev.A , 063630 (2010).[31] A. S. Bradley and B. P. Anderson, Phys. Rev. X ,041001 (2012).[32] M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S. Bradley, Phys. Rev. Lett. , 104501 (2013).[33] T. Bland, G. W. Stagg, L. Galantucci, A. W. Baggaley,and N. G. Parker, Phys. Rev. Lett. , 174501 (2018).[34] G. W. Stagg, N. G. Parker, and C. F. Barenghi, Phys.Rev. Lett. , 135301 (2017).[35] N. G. Berloff and C. Yin, J. Low Temp. Phys. , 187(2006).[36] H. Takeuchi, S. Ishino, and M. Tsubota, Phys. Rev. Lett. , 205301 (2010).[37] K. Fujimoto and M. Tsubota, Phys. Rev. A , 033642(2012); Phys. Rev. A , 053641 (2012); Phys. Rev. A ,063628 (2013); Phys. Rev. A , 013629 (2014); Phys.Rev. A , 033620 (2016).[38] M. Tsubota, Y. Aoki, and K. Fujimoto, Phys. Rev. A ,061601(R) (2013).[39] B. Villase˜nor, R. Zamora-Zamora, D. Bernal, and V.Romero-Roch´ın, Phys. Rev. A , 033611 (2014).[40] D. Kobyakov, A. Bezett, E. Lundh, M. Marklund, andV. Bychkov, Phys. Rev. A , 013631 (2014).[41] S. Kang, S. W. Seo, J. H. Kim, and Y. Shin, Phys. Rev.A , 053638 (2017).[42] M. Karl, B. Nowak, and T. Gasenzer, Sci. Rep. , 2394(2013); Phys. Rev. A , 063615 (2013).[43] K. Kudo and Y. Kawaguchi, Phys. Rev. A , 013630(2013); Phys. Rev. A , 053609 (2015).[44] J. Hofmann, S. S. Natu, and S. Das Sarma, Phys. Rev.Lett. , 095702 (2014).[45] S. De, D. L. Campbell, R. M. Price, A. Putra, B. M.Anderson, and I. B. Spielman, Phys. Rev. A , 033631(2014).[46] E. Nicklas, M. Karl, M. H¨ofer, A. Johnson, W. Mues-sel, H. Strobel, J. Tomkoviˇc, T. Gasenzer, and M. K.Oberthaler, Phys. Rev. Lett. , 245301 (2015).[47] L. A. Williamson and P. B. Blakie, Phys. Rev. Lett. ,025301 (2016); Phys. Rev. A , 023608 (2016); Phys.Rev. Lett. , 255301 (2017).[48] A. Bourges and P. B. Blakie, Phys. Rev. A , 023616(2017).[49] M. Pr¨ufer, P. Kunkel, H. Strobel, S. Lannig, D. Linne-mann, C.-M. Schmied, J. Berges, T. Gasenzer, and M.K. Oberthaler, Nature (London) , 217 (2018).[50] K. Fujimoto, R. Hamazaki, and M. Ueda, Phys. Rev.Lett. , 073002 (2018).[51] H. Takeuchi, Phys. Rev. A , 013617 (2018).[52] L. M. Symes, D. Baillie, and P. B. Blakie, Phys. Rev. A , 063618 (2018).[53] See, e.g., L. D. Landau and E. M. Lifshitz, Fluid Mechan-ics , 2nd ed. (Butterworth-Heinemann, Oxford, 1987),Chap. VII.[54] P. Ao and S. T. Chui, Phys. Rev. A , 4836 (1998).[55] R. A. Barankov, Phys. Rev. A , 013612 (2002).[56] B. Van Schaeybroeck, Phys. Rev. A , 023624 (2008);Phys. Rev. A , 065601 (2009).[57] S. Choi, S. A. Morgan, and K. Burnett, Phys. Rev. A ,4057 (1998).[58] See, e.g., C. J. Pethick and H. Smith, Bose-Einstein con-densation in dilute gases , 2nd ed., Chap. 12. (CambridgeUniv. Press, Cambridge, 2008).[59] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.Flannery,
Numerical Recipes , 3rd ed. (Cambridge Univ.Press, Cambridge, 2007).[60] See Supplemental Material athttp://link.aps.org/supplemental/ for movies of thedynamics. [61] M. R. Andrews, C. G. Townsend, H.-J. Miesner, D. S. Durfee, D. M. Kurn, and W. Ketterle, Sience275