Observation of bistable turbulence in quasi-two-dimensional superflow
OObservation of bistable turbulence in quasi-two-dimensional superflow
E. Varga, V. Vadakkumbatt, A. J. Shook, P. H. Kim, and J. P. Davis Department of Physics, University of Alberta, Edmonton, Alberta T6G 2E1, Canada
Turbulent flow restricted to two dimensions can spontaneously develop order on large scales,defying entropy expectations and in sharp contrast with turbulence in three dimensions wherenonlinear turbulent processes act to destroy large-scale order. In this work we report the observationof unusual turbulent behavior in steady-state flow of superfluid He—a liquid with vanishing viscosityand discrete vorticity—in a nearly two-dimensional channel. Surprisingly, for a range of experimentalparameters, turbulence is observed to exist in two bistable states. This bistability can be wellexplained by the appearance of large-scale regions of flow of opposite vorticity.
Chaotic motion of flowing fluids—turbulence—is oneof the most ubiquitous phenomena occurring in natureand is frequently encountered in everyday life. Typically,the turbulence that one encounters takes place in threedimensions (3D), however, two dimensional (2D) turbu-lence, while not perfectly realized in nature, is relevantto systems where motion in two dimensions dominatesover the third, such as large-scale flows in oceans [1], at-mospheres [2], soap bubbles [3] or liquid crystal films [4].The hallmark feature of turbulence in three dimensionsis the transfer of energy from large scales to small scalesin both classical [5] and quantum [6] fluids. This “Kol-mogorov cascade” can be understood as the splitting oflarge eddies in the flow into progressively smaller ones,until viscous damping dominates and dissipates the ki-netic energy of small scales into heat. Turbulence in 3Dtherefore acts to destroy any large scale ordering and,indeed, homogeneous and isotropic turbulence is an ex-cellent approximation in many cases.Restricting the flow of a classical fluid to 2D disruptsthis homogenizing behavior. Interestingly, the directionof the cascade of energy can be inversed [7, 8] and vortic-ity can coalesce into large eddies, thus spontaneously gen-erating large-scale order from forcing on smaller scales. Ifthe vorticity of the system is discrete (e.g., the quantizedvortices in Bose-Einstein condensates [9, 10] or the super-fluid phases of He and He [11], as opposed to continuousvorticity of classical fluids), one can treat the system asa ‘gas’ of point-like vortices, which can be analyzed us-ing the tools of statistical mechanics. In his pioneeringwork, Onsager [12] showed that such a gas can exist at ef-fectively negative temperatures, which would physicallymanifest as clusters of like-signed vortices (i.e., configu-rations with high energy and low entropy), similar to thelarge-scale eddies in 2D classical fluids.Sixty years after its prediction, the Onsager vortex gashas recently been observed: first using quantized vor-tices in Bose-Einstein condensates (BECs) [9, 10] andthen using a nanometer-thick film of superfluid He [11].These systems, however, contained only a small numberof vortices (
N <
50) and were allowed to decay freelyduring the experiment. Therefore open questions remainas to the robustness of this phenomenon in macroscopic systems with large number of vortices and in steady-state flows (regimes approached so far only in simula-tions [13, 14]). In this work, we study a forced andstrongly turbulent oscillatory flow in a µ m-thick slabof superfluid He with macroscopic (mm-scale) lateralsize. Turbulence in this system can exist in two nearly-degenerate bistable states, both different from the lami-nar (i.e., non-turbulent) state. The transitions betweenthese individual flow states are discontinuous, hysteretic,and a highly unusual “backward” transition from a less-turbulent to more-turbulent state upon decrease in veloc-ity is observed. We argue that these observations stemfrom quasi-2D physics and that both the bistability andbackward transition are naturally explained in terms ofspontaneous flow polarization, suggesting the presence oflarge-scale order.We study turbulence in superfluid He (He-II), whichbehaves as a mixture of two distinct fluid components[15]—an inviscid superfluid, where vorticity is restrictedto discrete quantized vortices, and a normal fluid, withcontinuous vorticity and finite viscosity. He-II has provento be a valuable testbed [16–19] for the study of turbu-lence with both continuous and discrete vorticity, as wellas the interactions between them.Here, oscillatory flow is excited inside a microfluidicHelmholtz resonator [20–22] immersed in He-II, whereflow is limited to nearly two dimensions by confinementin the vertical direction ( D = 1067 nm, Fig. 1(a)). A uni-form confinement, while somewhat larger than the thick-ness of previously used adsorbed films [11], avoids dissi-pative effects stemming from vortex-surface interactions[23]. Due to this strong confinement, only the superfluidcomponent of He-II can move (the normal componentbeing viscously clamped [21]). The resonator is micro-fabricated from single-crystal quartz and consists of acentral circular basin connected to a surrounding bath ofHe-II through two equal opposing channels of rectangu-lar cross-section (Fig. 1(a)). The capacitive driving andsensing of the flow (see [20–22] and [24] for more details)allows us to measure the relationship between the drivingpressure gradient and the fluid velocity in the channel ofthe resonator. The confinement used in this study was1067 nm, but qualitatively similar results were also ob- a r X i v : . [ c ond - m a t . o t h e r] J u l tained for 805 nm confinement (see [24]).A simulation of the fluid Helmholtz mechanical mode,Fig. 1(b), shows that the flow velocity is essentially con-fined to the two side channels. As the fluid flows intoor out of the channel, the sharp corners at the channelend induce a vorticity in the flow, as seen in Fig. 1(c,d).On the forward stroke (fluid flowing into the channel,Fig. 1(c)) vorticity is injected into the channel in a po-larized fashion. On the reverse stroke (Fig. 1(d)) vorticityis ejected and lost into the basin. The two side channelsare identical, thus any flow instabilities are likely to occurapproximately simultaneously.To study the dissipation in this quasi-2D flow we reso-nantly drive the Helmholtz mechanical mode and contin-uously increase or decrease the drive amplitude (no sig-nificant dependence on ramp rate was observed). Severalrepeated measurements of peak velocity as a function ofpeak applied pressure at nominally identical conditionsare shown in Figs. 2(a)-2(b). For small drives the behav-ior is linear (i.e., the flow is laminar). With increasingdrive, however, the measured velocity falls short of thevalue expected by extrapolating from the linear regime.That is, above a critical velocity the flow is damped by adrag with nonlinear dependence on velocity. The damp-ing in the linear regime is believed to be dominated by thethermoviscous effect [21], whereas the nonlinear dampingis predominantly due to presence of quantized vortices[37]. The transition to the nonlinear regime is hysteretic (a)(b) (c) w = l = D = 1067 nm -1 10normalized vorticity n o r m a li z e d r m s v e l o c i t y fl owdirection (d) FIG. 1. Polarized vorticity in a 2D Helmholtz resonator. (a)Sketch of the device. The central circular basin, which is onlyused to drive and sense the flow, is connected to the surround-ing bath of He-II through two side channels. (b) Simulationof the Helmholtz mechanical mode, showing the normalizedroot-mean-square velocity concentrated in the channels. (c,d)Vorticity is induced by sharp corners during the forward andreverse stroke of the mechanical mode; arrow indicates thesuperfluid flow direction. The positive and negative vorticesare well separated in space (cf. the vortex generation terms g ± in Eqs. 1,2) and the injection of vorticity into the channel isthus strongly polarized (cf. the vortex polarization generationterm g s in Eq. 4). c r i t i c a l v e l o c i t y ( m / s ) IIIIII v e l o c i t y ( m / s ) (b)(c) IIII II v e l o c i t y ( m / s ) (a) bistableturbulence FIG. 2. Bistable turbulence and critical velocities. (a) Themeasured flow velocity as a function of applied pressure gradi-ent for a range of temperatures. Darker curves show increas-ing pressure gradient, lighter decreasing (as indicated by thecurved arrows). (b) Detail of the laminar-to-turbulent tran-sition at 1.4 K. Blue curves correspond to increasing drive,orange to decreasing drive. Three critical velocities with dis-continuous jumps are apparent: “I” – transition from laminarto turbulent state upon increasing velocity; “II” – the un-usual backward transition into a more dissipative state upondecrease in velocity; “III” – transition from turbulent stateback to laminar flow. Above “I”, the flow can randomly tran-sition into the more dissipative state (not shown, see [24]). (c)The temperature dependence of the critical velocities of thethree types labeled in panel (b). Onset of the bistable tur-bulence coincides with the appearance of the critical velocity“II” of the backward transition. and is marked by a discontinuous jump in the velocity-pressure dependence. For temperatures below 1.7 K, thevelocity-pressure dependence in the nonlinear regime ran-domly follows one of two distinct and well-defined curves,i.e., the turbulence is bistable.The temperature dependence of the observed criticalvelocities (defined as the mean positions of the discon-tinuous jumps, see Fig. 2(b)) is shown in Fig. 2(c). Here,a new critical velocity—type “II”—appears below 1.7 K,which coincides with beginning of the bistable regime.In this bistable regime, as the flow velocity decreases theintermediate turbulent state with lower dissipation (i.e.,higher velocity at a given drive) destabilises but, ratherthan becoming laminar again, the flow transitions intothe state with stronger turbulence (i.e., lower velocity ata given drive), as can be seen in Fig. 2(b). This results ina highly unusual “backward” transition (critical velocity“II” in Fig. 2(b)) into a state with higher dissipation asthe flow velocity decreases .The microscopic confinement and large aspect ratio ofour flow suggests the use of a 2D theory such as theOnsager vortex gas model [7, 12]. However, our systemdeviates from the Onsager model in several importantaspects: it is dissipative, continuously driven, and theconfinement is large compared to the thickness of a quan-tized vortex ( ≈ − m). Since our measurements havebeen conducted at relatively high temperatures, mutualfriction will strongly attenuate any highly-curved vortexstructures. Therefore, for the sake of simplicity of themodelling, we assume that the majority of the vorticesin our system can be described by two populations withdefinite orientations, i.e., the vortices are approximatelypoint-like (see [24] for more detailed analysis). We note,however, that a population of vortices without definitepolarization (i.e., loops attached to a single wall) almostcertainly exists in our system. In quasi-2D modellingthese can be approximated as point-vortex dipoles. Fur-thermore, our experiment is sensitive to the total dissipa-tion, which is an integral quantity, and hence we cannotdirectly determine the presence of, e.g., negative vortextemperatures, for which we would need to know the po-sitions and signs of the vortices [38, 39]. However, thespatial separation of vortices of differing signs explainsthe observed bistability, hysteresis, and backward transi-tion.To show this, we construct a model for the num-ber of vortices in the system that captures the essentialphysics. A similar approach has been adapted for 2DBECs [40, 41] and 3D counterflow of He-II [42]. We modelthe time evolution (on time-scales long comparable tothe flow oscillation period) of positively- and negatively-oriented local vortex densities, n + and n − , respectivelyas ∂n + ∂t = an + + bn − − n + n − d + g + , (1) ∂n − ∂t = an − + bn + − n + n − d + g − . (2)Here, the terms on the right-hand-side correspond to re-moval of vortices by advection ( a < b > d > g ± > n = n + + n − and polarization s = ( n + − n − ) /n wehave ∂n∂t = ( a + b ) n − dn (1 − s ) + g, (3) ∂s∂t = − bs + 12 dns (1 − s ) + g s n , (4) v o r t e x nu m b e r p o l a r i z a t i o n n v o r t e x p o l a r i z a t i o n s (a) (b)(c) (d) s < 0 s > 0 g , gs ramp down FIG. 3. Bistability of the quasi-2D model, Eqs. 3, 4, for anexample set of parameters a = − b = 0 . d = 3, g = 2,and g s = 0 .
1. (a) Starting from initial conditions indicatedby the filled circles, the solution to equations 3, 4 approaches(indicated by the color of the trajectory) one of the stationarypoints shown as the blue stars. (b) Illustration of the two tur-bulent states. The large-scale polarization of the flow is eitheraligned ( s >
0) or anti-aligned ( s <
0) with the polarizationof the drive (cf. Fig. 1(c)). (c,d) The backward transitionin n (c) and s (d) during linear (in time) ramp-down of thegeneration parameters g and g s in the range 100 < t < s > n , however, s < s < s > where g = g + + g − and g s = [(1 − s ) g + − (1+ s ) g − ], whichwe take as our control parameters and assume their in-dependence of s and n (see [24] for details). The term g s represents the polarization of the drive, i.e., the sepa-ration of generation of positive and negative vortices onthe opposing corners of the channel (see Fig. 1). Thevortex densities n ± and vortex generation terms g ± arelocal, whereas only the total drag, determined by the to-tal number of vortices n , is measured in the experiment.As a first approximation we can replace the density n byits spatial average. The polarization s is anti-symmetricwith respect to the device axis (see Fig. 1) and its av-erage vanishes, assuming that the flow remains neutral.Therefore we decompose s ( r ) into a series of appropri-ate orthogonal modes s ( r ) = (cid:80) k s k ( r ). Truncating theexpansion after the leading term, we use Eq. 4 for calcu-lating the evolution of a single mode which captures thelarge-scale polarization of the vortex distribution.The dynamical system of Eqs. 3, 4 indeed has two sta-tionary solutions for certain choices of parameters thatdiffer in both s and n , as shown in Fig. 3(a). The essen-tial reason for the existence of two distinct steady states d nu m b e r o f v o r t i c e s n FIG. 4. The appearance of the bistable behavior controlledby d (other parameters remain unchanged). The parameter d is related to the cross-section for reconnection of collidingvortices. With decreasing temperature, the vortex motion isless damped thus making the vortices more curved. “Wig-gly” vortices occupy larger area and are therefore more likelyto reconnect. The appearance of the bistability with increas-ing d is therefore in qualitative agreement with the measuredtemperature dependence. of vortex number n is that g s (cid:54) = 0 (i.e., the drive is po-larized, see Fig. 1(c) and Fig. 3(b)), which lifts the de-generacy of s > s < g, g s linearly in time, we model the velocity ramp-down experiment. The result, shown in Figs. 3(c)- 3(d),reproduces the unusual backward transition observed ex-perimentally. The transition occurs due to destabilisa-tion of the s < n (i.e., the flowstate is robust enough to withstand the oppositely polar-ized drive). As the drive, and the overall vortex number,decreases, this state destabilizes and transitions into the s > n temporarily increases.Finally, the temperature dependence of the experimen-tal observations can be connected, for example, with theparameter d , which is related to the cross-section for re-connection of colliding vortices. This will increase withvortex deformation which, in turn, is expected to increasewith decreasing temperature [46]. As shown in Fig. 4, thebistability does indeed appear as d increases in qualita-tive agreement with the data.In conclusion, using a microfluidic Helmholtz resonatorwe have demonstrated a long-lived bistable turbulent be-havior in superfluid He restricted to quasi-2D channel,which exists below a certain critical temperature. In ad-dition, we observe an unusal backward transition wherethe flow transitions into a more dissipative state as theflow velocity decreases . The bistability, hysteresis andthe backward transition of the observed turbulence areunderstood in terms of a model of vortex density as an interplay between spontaneous flow ordering and polar-ization of turbulence generation. The proposed modelis, in principle, applicable to other systems with discretevorticity (e.g., BECs, superfluid He) if the generationof turbulence is in some manner polarized. An interest-ing question is whether similar behavior is possible incontinuous classical systems (indeed, random switchingbetween two degenerate flow configurations has been ob-served [47]). The backward transition is of particularinterest, as one usually expects turbulent fluctuations todecrease as the flow driving them is reduced. Consideringthat the driving mechanisms of, for example, atmosphericor oceanic flows—which are approximately 2D on largescales [2]—are typically not homogeneous and isotropic,the bistable behavior could have implications for weatherprediction, climate modelling, and atmospheres of gas gi-ants [48].This work was supported by the University of Alberta,Faculty of Science; the Natural Sciences and EngineeringResearch Council, Canada (Grants Nos. RGPIN-04523-16, DAS-492947- 16, and CREATE-495446-17); and theCanada Foundation for Innovation. We are grateful toG.G. Popowich for technical assistance and F. Souris forthe velocity calibration theory. [1] B. K. Arbic, K. L. Polzin, R. B. Scott, J. G. Richman,and J. F. Shriver,
J. of Phys. Oceanogr. , 283 (2013).[2] A. Pouquet and R. Marino, Phys. Rev. Lett. , 234501(2013).[3] M. Rivera, P. Vorobieff, and R. E. Ecke,
Phys. Rev. Lett. , 1417 (1998).[4] P. Tsai, Z. A. Daya, and S. W. Morris, Phys. Rev. Lett. , 084503 (2004).[5] U. Frisch. Turbulence: The legacy of A. N. Kolmogorov (Cambridge University Press, 1995).[6] N. Navon, C. Eigen, J. Zhang, R. Lopes, A. L. Gaunt, K.Fujimoto, M. Tsubota, R. P. Smith, and Z. Hadzibabic,
Science , 382 (2019).[7] R. H. Kraichnan and D. Montgomery,
Rep. Prog. Phys. , 547 (1980).[8] G. Boffetta and R. E. Ecke, Annu. Rev. Fluid Mech. ,427 (2012).[9] G. Gauthier, M. T. Reeves, X. Yu, A. S. Bradley, M. A.Baker, T. A. Bell, H. Rubinsztein-Dunlop, M. J. Davis,and T. W. Neely, Science , 1264 (2019).[10] S. P. Johnstone, A. J. Groszek, P. T. Starkey, C. J.Billington, T. P. Simula, and K. Helmerson,
Science ,1267 (2019).[11] Y. P. Sachkou, C. G. Baker, G. I. Harris, O. R. Stockdale,S. Forstner, M. T. Reeves, X. He, D. L. McAuslan, A. S.Bradley, M. J. Davis, and W. P. Bowen,
Science ,1480 (2019).[12] L. Onsager,
Nuovo Cimento , 279 (1949).[13] M. T. Reeves, T. P. Billam, X. Yu, and A. S. Bradley, Phys. Rev. Lett. , 184502 (2017).[14] M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S.Bradley,
Phys. Rev. Lett. , 104501 (2013). [15] D. R. Tilley and J. Tilley,
Superfluidity and Supercon-ductivity . (IoP, third edition, 1990).[16] C. F. Barenghi, L. Skrbek, and K. R Sreenivasan,
Proc.Natl. Acad. Sci. U.S.A , 4647 (2014)[17] P. Walmsley, D. Zmeev, F. Pakpour, and A. Golov,
Proc.Natl. Acad. Sci. U.S.A , 4691 (2014)[18] E. Fonda, K. R Sreenivasan, and D. P Lathrop,
Proc.Natl. Acad. Sci. U.S.A , 1924 (2019)[19] A. W. Baggaley, C. F. Barenghi, and Y. A. Sergeev,
Phys.Rev. E , 013002 (2014).[20] X. Rojas and J. P. Davis, Phys. Rev. B , 024503(2015).[21] F. Souris, X. Rojas, P. H. Kim, and J. P. Davis, Phys.Rev. Appl. , 044008 (2017).[22] A. J. Shook, V. Vadakkumbatt, P. Senarath Yapa,C. Doolin, R. Boyack, P. H. Kim, G. G. Popowich,F. Souris, H. Christani, J. Maciejko, and J. P. Davis, Phys. Rev. Lett. , 015301 (2020).[23] S. Forstner, Y. Sachkou, M. Woolley, G. I. Harris, X. He,W. P. Bowen, and C. G. Baker,
New J. Phys. , 053029(2019).[24] See appendix for details on experimental setup, pressuregradient and velocity calibration, discussion on two di-mensionality, further analysis and extension of the theo-retical model and additional data.[25] H. A. Kierstead, J. Low Temp. Phys. , 791 (1976).[26] D. J. Bishop and J. D. Reppy, Phys. Rev. Lett. , 1727(1978).[27] S. Jose Benavides and A. Alexakis, J. Fluid Mech. ,364 (2017).[28] S. Babuin, E. Varga, L. Skrbek, E. L´evˆeque, and P.-E.Roche,
EPL , 24006 (2014)[29] M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S.Bradley,
Phys. Rev. Lett. , 155302 (2015).[30] A. Pouquet, R. Marino, P. D. Mininni, and D. Rosenberg,
Phys. Fluids , 111108 (2017).[31] A. Duh, A. Suhel, B. D. Hauer, R. Saeedi, P. H. Kim,T. S. Biswas, J.P. Davis, J. Low Temp. Phys. , 31-39(2012).[32] K. W. Schwarz,
Phys. Rev. B , 5782-5804 (1985).[33] R. J. Donnelly, Quantized Vortices in Helium II . Cam-bridge University Press, 1991.[34] R. J. Donnelly and C. F. Barenghi,
J. Phys. Chem. Ref.Data. , , 1217 (1998).[35] G. W. Stagg, N. G. Parker, and C. F. Barenghi, Phys.Rev. Lett. , 135301 (2017).[36] E. J. Hopfinger, F. K. Browand, and Y. Gagne,
J. FluidMech. , 505 (1982).[37] J. Gao, W. Guo, S. Yui, M. Tsubota, and W. F. Vinen,
Phys. Rev. B , 184518 (2018).[38] A. J. Groszek, M. J. Davis, D. M. Paganin, K. Helmerson,and T. P. Simula, Phys. Rev. Lett. , 034504 (2018).[39] R. N. Valani, A. J. Groszek, and T. P. Simula,
New J.Phys. , 053038 (2018).[40] A. J. Groszek, T. P. Simula, D. M. Paganin, and K.Helmerson. Phys. Rev. A , 043614 (2016).[41] A. Cidrim, F. E.A. Dos Santos, L. Galantucci, V. S. Bag-nato, and C. F. Barenghi, Phys. Rev. A , 033651(2016).[42] E. Varga and L. Skrbek, Phys. Rev. B , 064507 (2018).[43] A. A. Berryman, Ecology , 1530 (1992).[44] A. J. Lotka, J. Phys. Chem. , 271 (1910).[45] H. Y. Shih, T. L. Hsieh, and N. Goldenfeld, Nat. Phys. , 245 (2016). [46] C. F. Barenghi, R. J. Donnelly, and W. F. Vinen, Phys.Fluids , 498 (1985)[47] M. L. Woyciekoski, L. A. M. Endres, A. V. de Paula andS. V.Mller Ocean Eng. , 106658 (2020)[48] R. M. B. Young and P. L. Read,
Nat. Phys. , 1135(2017). SUPPLEMENTARY INFORMATIONEXPERIMENTAL SETUP AND SUPERFLUID VELOCITY CALCULATION.
Flow in the Helmholtz resonators is driven and sensed capacitively using two aluminum electrodes deposited on thetop and bottom wall of the device, forming a parallel plate capacitor (see [21] for details on the fabrication process).An alternating voltage of amplitude U applied to the electrodes of the device causes a periodic deformation of thewalls of the basin due to electrostatic force, which pushes the superfluid in and out of the basin through the two sidechannels and into the bath, thus driving the Helmholtz mode. The Helmholtz resonance is observed as a periodicvariation of the capacitance of the device.The resonator is wired in a bridge circuit shown in Fig. SM1, balanced to the capacitance C of the resonator at rest.Change in this capacitance caused by the Helmholtz resonance results in bridge imbalance and a current I throughthe detector G. An example of the resulting spectrum of two resonators wired in parallel is shown in Fig. SM2.The current is first amplified by a transimpedance amplifier (Stanford Research SR570) and then measured by a(Zurich Instruments HF2LI) lock-in amplifier referenced to the frequency of excitation U . A standard 9V battery isused as the source of the bias voltage ( U B = 9 . G U B C I C b U Basin ChannelAl electrodeBTBT
FIG. SM1. Measurement scheme. Aluminum electrodes of the Helmholtz resonator form a parallel plate capacitor of capacitance C biased by the battery U B . The resonator is wired as one arm of a capacitance bridge (the other arm being the capacitor C b ) which is balanced such that current I through the detector G is approximately zero when flow of the helium is negligible.When the Helmholtz mechanical mode is excited by the oscillating voltage U , the oscillating pressure in the basin changes thecapacitance of the device and thus a nonzero current through the detector G. The bridge circuit is isolated from the batteryvoltage U B by the two bias tees (BT). The transformer ratio is 1:1 on the resonator arm and adjustable for the C b arm. In this section we first analyze the equation of motion of the fluid in channels, approximated as a mass on a spring,and derive the relationship between the oscillating driving voltage U and the pressure gradient in the channels. Next,we calculate the superfluid velocity in the channels from the current I through the detector. Helmholtz equation of motion. –
We derive the equation of motion of the superfluid in the Helmholtz resonatorwith explicit drive and damping forces. We approximate the flow in the channel as a mass on a spring, which isdisplaced by distance y (positive in the direction away from the basin). The average displacement of the plates of thebasin is denoted by x (positive when the basin contracts). We begin by calculating the change in total density of thefluid inside the basin as a response to the mean deformation of the basin x and the displacement of the superfluidinside the channel y : δρ = δ (cid:18) M B V B (cid:19) = δM B V B − M B V B δV B = 1 V B ( − aρ s y + 2 ρAx ) . (5)Here M B = ρV B is the mass of the fluid inside the basin, V B = AD is the volume of the basin ( A being its area and D the confinement) and from the assumption that only the superfluid moves δM B = − aρ s y ( a = wD is the cross-sectional area of the channel; the factor of 2 comes from the two channels), and δV B = − Ax is the change in basinvolume due to motion of the plates. A change in density corresponds to a change in pressure via the compressibility χ , δρ = ρχδP , or δP = 1 ρχV B (2 ρAx − aρ s y ) . (6)Balancing forces on the plate (neglecting its inertia) yields F es = 12 k p x + AδP, (7)
800 1000 1200 1400 1600frequency (Hz)0.020.030.040.050.060.070.080.090.10 n o r m a li z e d d e t e c t o r c u rr e n t ( A / V ) d r i v e v o l t a g e ( V ) FIG. SM2. The spectrum of two Helmholtz resonators wired in parallel measured at 1.4 K and normalized by the drive voltageamplitude U . The fact that the normalized peaks do not overlap for increasing drives indicates nonlinear dissipation. Thehigher-frequency peak corresponds to the device with D = 1067 nm and the lower-frequency one to D = 805 nm (see Eq. 12). where F es = C U / (2 D ) is the electrostatic force between the parallel plates of the capacitor formed by the circularelectrodes in the basin, U = U B + U is the total applied voltage and k p = 2 . × N/m [21] is the stiffness of thesubstrate deflection (note that this is double that in Ref. [21], where the stiffness refers to deflection of both plates inparallel). Expressing x from Eq. 7 and substituting back in to Eq. 6 yields δP = k p ρ ( χV B k p + 4 A ) (cid:20) ρAk p F es − aρ s y (cid:21) . (8)The superfluid inside the channel is accelerated by the pressure ρ s al ¨ y = ρ s ρ aδP − F f , (9) ρ s al ¨ y = ρ s k p aρ ( χV B k p + 4 A ) (cid:20) ρAk p F es − aρ s y (cid:21) − F f , (10)where we included a friction force F f = alρ s ζ ˙ y . Here ζ is a friction parameter with units of frequency but will remainotherwise unspecified for now. Rearranging,¨ y + 2 ρ s k p alρ ( χV B k p + 4 A ) y + ζ ˙ y = 1 ρ Al ( χV B k p + 4 A ) F es , (11)from which the resonance frequency follows ω = 2 alρ ρ s ρ k p A (1 + Σ) , (12)where Σ = χDk p / (4 A ).Finally, the driving pressure gradient, the quantity shown on the x-axis in Fig. 2A,B of the main text, is given by δPl = 4 A ( χV B k p + 4 A ) l F es = 4 AC U B ( χV B k p + 4 A ) Dl U , (13)where we take only the component of the force F es on resonance with the Helmholtz mode ( U being the AC drive), F reses = C U U B /D . Calculation of velocity from detector current. –
Whenever the capacitance bridge in the measurement circuitshown in Fig. SM1 becomes imbalanced, current will flow through the detector. Assuming that the bridge is tunedto the total capacitance of the devices at rest, the current through the detector at the frequency of the drive (theonly component detected by the lock-in amplifier) is given only by the oscillation of the device capacitance due to theHelmholtz resonance, I = d( CU )d t = U B d C d t = U B d C d y ˙ y, (14)where U B is the bias voltage. The change of capacitance with superfluid displacement in the channel can be writtenas d C d y = C ε d ε d ρ d ρ d y + 2 C D d x d y , (15)where C = ε εA el /D is the capacitance of an undisturbed device with ε , ε being the vacuum permittivity anddielectric constant of helium, respectively, and A el the area of the electrodes. The spacing between electrodes is givenby h = D − x , hence the second term in Eq. 15.Neglecting the dependence of polarizability of helium on density [25] and using the Clausius-Mossoti relation wecan estimate the change in dielectric constant as d ε d ρ = ε − ρ . (16)The change of density with superfluid displacement can be calculated directly from Eq. 8 and using d ρ = ρχ d P ,d ρ d y = − aρ s V B . (17)where Σ = χDk p / (4 A ). Differentiating Eq. 5 with respect to y and putting the result equal to Eq. 17 results in anequation for d x/ d y and yields d x d y = aρ s Aρ (cid:18) − (cid:19) . (18)Inserting Eq. 18, Eq. 17 and Eq. 16 back into Eq. 15 yields g ≡ C d C d y = 2 aρ s V B ρ (cid:18) − ε − ε Σ (cid:19)
11 + 2Σ . (19)Finally, the flow velocity is calculated as ˙ y = IgV B C , (20)assuming that the background has been subtracted from I . Two-dimensionality of turbulence
To what extent can the studied flow be considered 2D? The thickness of the flow channel is too large (i.e., D is muchlarger than the coherence length, ≈ ˚A) for finite-size effects of 2D superfluidity to be relevant [26]. 2D turbulence,on the other hand, requires only that the fluctuating velocity is restricted to 2D. This is essentially controlled by thechannel aspect ratio and damping of the self-induced vortex motion.Turbulence in He-II, especially when forced by a pressure gradient in large systems, typically behaves quasi-classically—as a classical liquid with effective viscosity [28]. Thus we first verify that the turbulence could be consideredtwo-dimensional based on classical fluid dynamics criteria.The forced superflow induced by the Helmholtz resonance is naturally 2D, however, the device geometry (specifically,the sharp corners near where the basin connects to the channel) induces shear on the scale of the channel width w = 1 . w/D offorcing to confinement scale. Specifically, for w/D (cid:38) √ Re the turbulence develops the 2D inverse energy cascade.Here, Re is the Reynolds number, which we define for our system using the effective quasi-classical viscosity He-II[28] ν eff ≈ . κ as Re = wv s /ν eff . In our experiments w/D ≈ D = 1067 nm and the highest experimentallyachieved √ Re ≈ a ≈ . v i ( R ) = κ πR (cid:20) log (cid:18) Ra (cid:19) − (cid:21) , (21)and, for stationary normal fluid, the change in radius is given by [33]˙ R = α [ V s ( t ) − v i ( R )] , (22)where α is the mutual friction constant [34] and V s ( t ) = V s0 sin Ω t is the imposed superflow. We numerically integratethe evolution of R for a range of initial radii R = R ( t = 0) and velocity amplitudes, terminating the calculationwhen either R ≈ R ≈ D = 1 µ m and the loop reconnects with the opposingwall, thus transforming into a vortex dipole. As shown in Fig. SM3(b), the typical lifetime t ∗ of half-rings for theparameters typical of our experiment is shorter than the flow oscillation period T , reaching, at most, about 0.6 T for a very specific choice of parameters. Vortex loops attached to a surface are thus short-lived transient objects.Creation and expansion of these loops is a likely scenario for vortex splitting and unpolarized injection, which featurein the quasi-2D model of Eqs. (3,4) of the main text, discussed further in the next section.Note, however, that we neglected the effects of the opposing wall on the self-induced velocity of the ring. Thiswill cause the vortex to deform and be attracted to the opposing wall, thus slightly altering the lifetime. Changingthe phase of the oscillating flow either does not significantly influence the outcome or causes loops of all sizes toquickly decay. Changing the angle between the plane of the loop and flow velocity would result in a somewhat morecomplicated transient flow, which is, however, unlikely to terminate in a significantly different manner.The calculation outlined above is not valid for vortex loop radii R comparable with the surface roughness b ≈ al. [35] studied a flow close to an irregular surface using a simulation of vortices in a Bose-Einstein condensate in thezero-temperature limit using the Gross-Pitaevskii equation (GPE). In that work, a dense layer of vortices was foundin the rough landscape of the surface sustained by intrinsic nucleation of vortices on the protruding peaks of thesurface. It is possible that such a dense boundary layer exists in our case as well, however, results obtained using GPEought to be adopted with caution for helium at finite temperatures. Intrinsic nucleation of vortices in He-II requiressignificantly higher velocities and mutual friction at finite temperatures will strongly damp any small, highly-curvedvortex structures. Regardless, this boundary layer is expected to be confined to within the scale of surface roughness[35] which in our case is about 0.1% of the confinement thus making it unlikely to be a significant contribution to theobserved macroscopic drag.Finally, the vortices connecting the two confining walls can, in principle, deform arbitrarily on the scale of theconfinement, D ≈ a . We estimate the dynamical importance of these deformations by comparing their typicalrate of decay to the time scale of their forcing, i.e., the flow oscillation period. Assume that the vertical modes offlow, mediated by the vortex deformation, take the form of a cascade of Kelvin waves—helical wave modes on vortices[33]. The decay rate of a Kelvin wave mode of wave vector k is τ = κ/ παk [46]. The smallest admissible k ≈ π/D results for T = 1 . D = 1067 nm in τ ≈ µ s. Increasing temperature will decrease τ . The decay rate ofKelvin waves is thus significantly faster than the time scale of their pumping (i.e., flow period, which is of the orderof 1 ms) and comparable to the inverse frequency of the Kelvin mode itself [33], i.e., no Kelvin wave cascade is likelyto develop along the individual vortices since the largest scales are already in the dissipative range. Other modes of0 v e l o c i t y a m p li t u d e ( m / s ) t * / T (a) (b) v FIG. SM3. (a) Configuration of the flow and vortex for calculation of the lifetime of the loop. Under the imposed oscillatoryflow, the loop can either annihilate or expand, reconnect with the opposite wall and thus create a dipole pair of vorticesspanning the confinement. (b) Time, t ∗ relative to the flow oscillation period T , for a half ring to grow to the size of thechannel (and thus reconnect with opposite wall to form a vortex dipole, red scale) or to completely annihilate (blue scale) inoscillating superfluid flow of varying amplitudes. Flow frequency was assumed to be f = 1200 Hz ( T = 1 /f = 0 .
83 ms) andmutual friction constant α = 0 .
034 (corresponding to approximately 1.3 K). Higher α (temperature) results in typically shorterlifetime. vortex deformation (e.g., solitons [36]), which cannot be decomposed to Kelvin waves, are possible. However, sincethe local velocity of the deformed line, and thus its decay rate mediated by mutual friction, are primarily determinedby the local curvature, we expect the decay of these deformations to be comparable to that of the Kelvin waves. Theamplitude of thermally excited Kelvin waves is also expected to be negligible [46]. We therefore consider the vorticesin our system which span the confinement to be point-like. Decreasing the temperature, particularly below 1 K, wouldsuppress the mutual friction damping and allow the vortices to deform strongly. Therefore we expect the turbulenceto cease to be 2D-like at sufficiently low temperatures. Vortex density model
Discrete quantized vortices are transported by flow similar to how vorticity is transported in classical 2D flow [7]: ∂n ± ∂t + ( v s · ∇ ) n ± = ( n + + n − ) b ( v s ) − dn + n − + g ± , (23)where the terms on the right hand side correspond, respectively, to splitting, decay by collision and generation ofvortices by the external drive. Note that this is not simply a passive scalar transport with source terms, since v s depends on the vortex distribution. We expect the splitting rate b to depend on velocity, possibly exhibiting criticalbehavior itself. For simplicity, however, we take b to be constant since in a high-velocity regime it is likely to bedominated by the flow oscillation period, which is independent of velocity.Averaged over the flow oscillation period, the advection term ( v s · ∇ ) n ± will have no effect on the vortex density farfrom the system boundaries. In the region near the boundaries, however, some vortices will be transported toward thewall and annihilated, i.e., the average effect of the advection is to reduce the vortex number. Since the vortex densitywill vary on the scale of the channel width, we approximate the gradient term as ( v s · ∇ ) n ± ≈ v s n ± /w . Putting a = b − pv s /w , where p characterises the inhomogeneous distribution of vortices throughout the channel, we recoverEqs. 1,2 of the main text. The assumption of velocity-independent b and d limits the applicability of the model toturbulent states at relatively high velocity and makes it unsuitable for modelling transition to turbulence from thelaminar state or predicting the scaling of vortex number with velocity.Following from Eqs. 1, 2 of the main text, the total vortex density n = n + + n − and polarization s = ( n + − n − ) /n obey d n d t = ( a + b ) n − dn (1 − s ) + g, (24)1 W /2- W /2- L /2 L /2 uv FIG. SM4. Simplified geometry for the decomposition of s into orthogonal modes. and d s d t = − bs + 12 dns (1 − s ) + g s n , (25)where we grouped all terms depending on g ± to new terms g = g + + g − and g s = (1 − s ) g + − (1 + s ) g − . Thegeneration terms g ± in Eqs. 1, 2 of the main text represent extrinsic or intrinsic nucleation of vortices and are likelyto be concentrated near the sharp corners connecting the basin and the channel. The vortices generated at theseedges in a polarized configuration are advected into the channel where they contribute to the observed drag. Near thecorners, however, the polarization s will likely be dominated by the instantaneous flow and, averaged over the flowoscillation period, s ≈ g s ≈ g + − g − , independent of s . For simplicity we adopt g and g s as independentcontrol parameters, rather than g ± . It should be noted, however, that it is the assumption of s -independent g s thatallows for bistable solutions.The equations above are assumed to be local, but spatially averaged quantities are required for comparison withthe experiment. The total vortex density n is an always positive quantity and thus, to a first approximation, can bereplaced by its spatial average. The vortex polarization s , on the other hand, has a vanishing average since we assumethat the flow will remain on average neutral.To connect Eq. 25 to averaged quantities, let us consider the simplified device geometry shown in Fig. SM4. Thebasin is removed and a single channel runs through the entire length of the device, but otherwise we assume generalflow features similar to the real device (e.g., flow direction, behavior of the generation terms g, g s ). From the symmetryof the problem, s is anti-symmetric with respect to mirroring about either of the axes and thus can be decomposedinto orthogonal modes as s ( u, v ) = (cid:88) k,l s kl sin (cid:18) πkL u (cid:19) sin (cid:18) πlW v (cid:19) , (26)and the spatial average is then given by (cid:104) s (cid:105) = 1 / (cid:80) kl | s kl | . In the actual device geometry the modes in theexpansion Eq. 26 will be more complicated but could, in principle, be constructed by a suitable transformation of therectangular domain of Fig. SM4 onto the actual device geometry. Truncating the expansion at the lowest s mode(which is likely to be the dominant term in the generation g s ) allows us to essentially use Eqs. 24, 25 as they are andrecover the results from the main text.In principle higher modes s kl can be considered, where Eq. 25 would be replaced by a set of equations for each modecoupled through nonlinear terms. The generation term g s is unlikely to have a single-mode decomposition and thenonlinear terms (in s ) in Eq. 25 will excite higher modes at the expense of lower modes. This picture is fully consistentwith the forward enstrophy (quadratic integral of vorticity) cascade of classical 2D turbulence [7]. Higher modes willagain exhibit near-degeneracy of the s kl > s kl < g s and,possibly, by the lower-lying modes through the nonlinear terms. This will result in more general multi-stability of themean vortex number n , as illustrated in Fig. SM5b. As the flow velocity or drive increases, the system will randomlyselect either the s > s <
50 55 60 65 70 75 80pressure gradient (kPa/m)0.800.850.900.951.00 v e l o c i t y ( m / s ) n =0 s >0 s <0 s >0 s <0 s >0 s <0laminar turbulent (a) (b) ... increasing drive FIG. SM5. (a) Three distinct turbulent branches of the pressure-velocity curve at 1.4 K (the laminar regime, see Fig. 2a,b ofthe main text, is not shown in this plot). The multi-stable behaviour hints at the involvement of higher modes of the vortexpolarization. (b) An illustration of a tree of multi-stable states generated by the Eqs. 24, 25. v e l o c i t y ( m / s ) FIG. SM6. Velocity-pressure gradient relationship for the 805 nm confinement, under identical conditions as 1067 nm confine-ment shown in Fig. 2 of the main text and Fig. SM7. Darker curves show increasing pressure gradient, lighter decreasing. Thebistability is observed here in a weaker form between 1.6 and 1.8 K.
Comparison of turbulence in 805 nm and 1067 nm confinements
The velocity-pressure gradient curves for 805 nm confinement (shown in Fig. SM6) were measured in parallel withthe 1067 nm confinement (shown in Fig. SM7 and Fig. 2 of the main text) under identical conditions. The bistabilityis again present in the 805 nm confinement, but in a weaker form and in the temperature range of 1.6–1.8 K. Thebistability also tends to be suppressed at high velocities.Within the model of Sec. the bistability can be destroyed in several ways, apart from the already discussed3 v e l o c i t y ( m / s ) FIG. SM7. Velocity-pressure gradient relationship for the 1067 nm, same data as in Fig. 2 of the main text shown for a widerrange of the pressure gradient. Darker curves show increasing pressure gradient, lighter decreasing. A few random transitionsfrom the less-dissipative to the more-dissipative state are visible for the bistable turbulence at
T < . temperature dependence of the decay parameter d . For example, increasing the splitting rate b above a critical valuewill result in a single solution with s ≈ g s above a critical value will destabilise the solution with sign of s opposite to the sign of g s (i.e., opposing polarizationwill be overwhelmed by the strong drive). Additionally, the confinement of the device likely affects the d parameteras well – the smaller 805 nm device allows for lesser lateral deformation of the vortices and hence lowers the effectivecross-section for collision, thus reducing d , which tends to reduce the bistability.In fact the bistability is not necessarily destroyed completely. If the environmental disturbances (e.g., vibrationof the cryostat) are non-negligible compared to the relative stability of the less-stable state (controlled, for example,by the d parameter), the flow would stochastically transition to the more stable state whenever a sufficiently strongfluctuation randomly occurs. This scenario is consistent with the fact that the transition between the two turbulentstates in the temperature range 1.6–1.8 K for the 805 nm confinement does not appear to have a well defined criticalvelocity.The device-dependence of the dd