# Observation of bistable turbulence in quasi-two-dimensional superflow

###### Abstract

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 where nonlinear turbulent processes act to destroy large-scale order. In this work we report the observation of unusual turbulent behavior in steady-state flow of superfluid He—a liquid with vanishing viscosity and discrete vorticity—in a nearly two-dimensional channel. Surprisingly, for a range of experimental parameters, turbulence is observed to exist in two bistable states. This bistability can be well explained by the appearance of large-scale regions of flow of opposite vorticity.

Chaotic motion of flowing fluids—turbulence—is one of the most ubiquitous phenomena occurring in nature and is frequently encountered in everyday life. Typically, the turbulence that one encounters takes place in three dimensions (3D), however, two dimensional (2D) turbulence, while not perfectly realized in nature, is relevant to systems where motion in two dimensions dominates over the third, such as large-scale flows in oceans Arbic2013 , atmospheres Pouquet2013 , soap bubbles Rivera1998 or liquid crystal films Tsai2004 . The hallmark feature of turbulence in three dimensions is the transfer of energy from large scales to small scales in both classical Frisch_book and quantum Navon2019 fluids. This “Kolmogorov cascade” can be understood as the splitting of large eddies in the flow into progressively smaller ones, until viscous damping dominates and dissipates the kinetic energy of small scales into heat. Turbulence in 3D therefore acts to destroy any large scale ordering and, indeed, homogeneous and isotropic turbulence is an excellent approximation in many cases.

Restricting the flow of a classical fluid to 2D disrupts this homogenizing behavior. Interestingly, the direction of the cascade of energy can be inversed Kraichnan1980 ; Boffetta2012 and vorticity can coalesce into large eddies, thus spontaneously generating large-scale order from forcing on smaller scales. If the vorticity of the system is discrete (e.g., the quantized vortices in Bose-Einstein condensates Gauthier2019 ; Johnstone2019 or the superfluid phases of He and He Sachkou2019 , as opposed to continuous vorticity of classical fluids), one can treat the system as a ‘gas’ of point-like vortices, which can be analyzed using the tools of statistical mechanics. In his pioneering work, Onsager Onsager1949 showed that such a gas can exist at effectively negative temperatures, which would physically manifest as clusters of like-signed vortices (i.e., configurations with high energy and low entropy), similar to the large-scale eddies in 2D classical fluids.

Sixty years after its prediction, the Onsager vortex gas has recently been observed: first using quantized vortices in Bose-Einstein condensates (BECs) Gauthier2019 ; Johnstone2019 and then using a nanometer-thick film of superfluid He Sachkou2019 . These systems, however, contained only a small number of vortices () and were allowed to decay freely during the experiment. Therefore open questions remain as 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 simulations Reeves2017 ; Reeves2013a ). In this work, we study a forced and strongly turbulent oscillatory flow in a µm-thick slab of superfluid He with macroscopic (mm-scale) lateral size. Turbulence in this system can exist in two nearly-degenerate bistable states, both different from the laminar (i.e., non-turbulent) state. The transitions between these individual flow states are discontinuous, hysteretic, and a highly unusual “backward” transition from a less-turbulent to more-turbulent state upon decrease in velocity is observed. We argue that these observations stem from quasi-2D physics and that both the bistability and backward transition are naturally explained in terms of spontaneous flow polarization, suggesting the presence of large-scale order.

We study turbulence in superfluid He (He-II), which behaves as a mixture of two distinct fluid components Tilley_book —an inviscid superfluid, where vorticity is restricted to discrete quantized vortices, and a normal fluid, with continuous vorticity and finite viscosity. He-II has proven to be a valuable testbed Barenghi2014c ; Walmsley2014a ; Fonda2019 ; Baggaley2014b for the study of turbulence with both continuous and discrete vorticity, as well as the interactions between them.

Here, oscillatory flow is excited inside a microfluidic Helmholtz resonator Rojas2015 ; Souris2017 ; Shook2020 immersed in He-II, where flow is limited to nearly two dimensions by confinement in the vertical direction ( nm, Fig. 1(a)). A uniform confinement, while somewhat larger than the thickness of previously used adsorbed films Sachkou2019 , avoids dissipative effects stemming from vortex-surface interactions Forstner2019 . Due to this strong confinement, only the superfluid component of He-II can move (the normal component being viscously clamped Souris2017 ). The resonator is microfabricated from single-crystal quartz and consists of a central circular basin connected to a surrounding bath of He-II through two equal opposing channels of rectangular cross-section (Fig. 1(a)). The capacitive driving and sensing of the flow (see Rojas2015 ; Souris2017 ; Shook2020 and SM for more details) allows us to measure the relationship between the driving pressure gradient and the fluid velocity in the channel of the resonator. The confinement used in this study was 1067 nm, but qualitatively similar results were also obtained for 805 nm confinement (see SM ).

A simulation of the fluid Helmholtz mechanical mode, Fig. 1(b), shows that the flow velocity is essentially confined to the two side channels. As the fluid flows into or out of the channel, the sharp corners at the channel end 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 polarized fashion. On the reverse stroke (Fig. 1(d)) vorticity is ejected and lost into the basin. The two side channels are identical, thus any flow instabilities are likely to occur approximately simultaneously.

To study the dissipation in this quasi-2D flow we resonantly drive the Helmholtz mechanical mode and continuously increase or decrease the drive amplitude (no significant dependence on ramp rate was observed). Several repeated measurements of peak velocity as a function of peak applied pressure at nominally identical conditions are shown in Figs. 2(a)-2(b). For small drives the behavior is linear (i.e., the flow is laminar). With increasing drive, however, the measured velocity falls short of the value expected by extrapolating from the linear regime. That is, above a critical velocity the flow is damped by a drag with nonlinear dependence on velocity. The damping in the linear regime is believed to be dominated by the thermoviscous effect Souris2017 , whereas the nonlinear damping is predominantly due to presence of quantized vortices Gao2018a . The transition to the nonlinear regime is hysteretic and is marked by a discontinuous jump in the velocity-pressure dependence. For temperatures below 1.7 K, the velocity-pressure dependence in the nonlinear regime randomly follows one of two distinct and well-defined curves, i.e., the turbulence is bistable.

The temperature dependence of the observed critical velocities (defined as the mean positions of the discontinuous 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 the intermediate turbulent state with lower dissipation (i.e., higher velocity at a given drive) destabilises but, rather than becoming laminar again, the flow transitions into the state with stronger turbulence (i.e., lower velocity at a given drive), as can be seen in Fig. 2(b). This results in a highly unusual “backward” transition (critical velocity “II” in Fig. 2(b)) into a state with *higher* dissipation as the flow velocity *decreases*.

The microscopic confinement and large aspect ratio of our flow suggests the use of a 2D theory such as the Onsager vortex gas model Onsager1949 ; Kraichnan1980 . However, our system deviates from the Onsager model in several important aspects: it is dissipative, continuously driven, and the confinement is large compared to the thickness of a quantized vortex ( m). Since our measurements have been conducted at relatively high temperatures, mutual friction will strongly attenuate any highly-curved vortex structures. Therefore, for the sake of simplicity of the modelling, we assume that the majority of the vortices in our system can be described by two populations with definite orientations, i.e., the vortices are approximately point-like (see SM for more detailed analysis). We note, however, that a population of vortices without definite polarization (i.e., loops attached to a single wall) almost certainly exists in our system. In quasi-2D modelling these can be approximated as point-vortex dipoles. Furthermore, our experiment is sensitive to the total dissipation, which is an integral quantity, and hence we cannot directly determine the presence of, e.g., negative vortex temperatures, for which we would need to know the positions and signs of the vortices Groszek2018a ; Valani2018 . However, the spatial separation of vortices of differing signs explains the observed bistability, hysteresis, and backward transition.

To show this, we construct a model for the number of vortices in the system that captures the essential physics. A similar approach has been adapted for 2D BECs Groszek2016 ; Cidrim2016 and 3D counterflow of He-II Varga2018 . We model the time evolution (on time-scales long comparable to the flow oscillation period) of positively- and negatively-oriented local vortex densities, and , respectively as

(1) |

(2) |

Here, the terms on the right-hand-side correspond to removal of vortices by advection (), creation of new vortices by splitting of seed vortices (), annihilation of a pair of vortices of opposite orientation (), and creation of vortices by large scale shear (). We note in passing that these equations are similar to the Lotka-Volterra, or predator-prey, equations used to model population dynamics in ecology Berryman1992 , oscillatory chemical reactions Lotka1910 or, indeed, the transition to turbulence Shih2016 . Restricting the model to total vortex density and polarization we have

(3) |

(4) |

where and , which we take as our control parameters and assume their independence of and (see SM for details). The term represents the polarization of the drive, i.e., the separation of generation of positive and negative vortices on the opposing corners of the channel (see Fig. 1). The vortex densities and vortex generation terms are local, whereas only the total drag, determined by the total number of vortices , is measured in the experiment. As a first approximation we can replace the density by its spatial average. The polarization is anti-symmetric with respect to the device axis (see Fig. 1) and its average vanishes, assuming that the flow remains neutral. Therefore we decompose into a series of appropriate orthogonal modes . Truncating the expansion after the leading term, we use Eq. 4 for calculating the evolution of a single mode which captures the large-scale polarization of the vortex distribution.

The dynamical system of Eqs. 3, 4 indeed has two stationary solutions for certain choices of parameters that differ in both and , as shown in Fig. 3(a). The essential reason for the existence of two distinct steady states of vortex number is that (i.e., the drive is polarized, see Fig. 1(c) and Fig. 3(b)), which lifts the degeneracy of and stationary solutions, when they exist.

By starting the evolution of the system from either of the stationary solutions and switching off the generation terms 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 experimentally. The transition occurs due to destabilisation of the polarization state, which is stabilised only at sufficiently high vortex densities (i.e., the flow state is robust enough to withstand the oppositely polarized drive). As the drive, and the overall vortex number, decreases, this state destabilizes and transitions into the state. When the flow and drive polarizations are aligned, fewer vortices are annihilated and thus the vortex density temporarily increases.

Finally, the temperature dependence of the experimental observations can be connected, for example, with the parameter , which is related to the cross-section for reconnection of colliding vortices. This will increase with vortex deformation which, in turn, is expected to increase with decreasing temperature Barenghi1985a . As shown in Fig. 4, the bistability does indeed appear as increases in qualitative agreement with the data.

In conclusion, using a microfluidic Helmholtz resonator we have demonstrated a long-lived bistable turbulent behavior in superfluid He restricted to quasi-2D channel, which exists below a certain critical temperature. In addition, we observe an unusal backward transition where the flow transitions into a *more* dissipative state as the flow velocity *decreases*. The bistability, hysteresis and the backward transition of the observed turbulence are understood in terms of a model of vortex density as an interplay between spontaneous flow ordering and polarization of turbulence generation. The proposed model is, in principle, applicable to other systems with discrete vorticity (e.g., BECs, superfluid He) if the generation of turbulence is in some manner polarized. An interesting question is whether similar behavior is possible in continuous classical systems (indeed, random switching between two degenerate flow configurations has been observed Woyciekoski2020 ). The backward transition is of particular interest, as one usually expects turbulent fluctuations to decrease as the flow driving them is reduced. Considering that the driving mechanisms of, for example, atmospheric or oceanic flows—which are approximately 2D on large scales Pouquet2013 —are typically not homogeneous and isotropic, the bistable behavior could have implications for weather prediction, climate modelling, and atmospheres of gas giants Young2017 .

This work was supported by the University of Alberta, Faculty of Science; the Natural Sciences and Engineering Research Council, Canada (Grants Nos. RGPIN-04523-16, DAS-492947- 16, and CREATE-495446-17); and the Canada Foundation for Innovation. We are grateful to G.G. Popowich for technical assistance and F. Souris for the velocity calibration theory.

## References

- (1) B. K. Arbic, K. L. Polzin, R. B. Scott, J. G. Richman, and J. F. Shriver, J. of Phys. Oceanogr. 43, 283 (2013).
- (2) A. Pouquet and R. Marino, Phys. Rev. Lett. 111, 234501 (2013).
- (3) M. Rivera, P. Vorobieff, and R. E. Ecke, Phys. Rev. Lett. 81, 1417 (1998).
- (4) P. Tsai, Z. A. Daya, and S. W. Morris, Phys. Rev. Lett. 92, 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 366, 382 (2019).
- (7) R. H. Kraichnan and D. Montgomery, Rep. Prog. Phys. 43, 547 (1980).
- (8) G. Boffetta and R. E. Ecke, Annu. Rev. Fluid Mech. 44, 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 364, 1264 (2019).
- (10) S. P. Johnstone, A. J. Groszek, P. T. Starkey, C. J. Billington, T. P. Simula, and K. Helmerson, Science 364, 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 366, 1480 (2019).
- (12) L. Onsager, Nuovo Cimento 6, 279 (1949).
- (13) M. T. Reeves, T. P. Billam, X. Yu, and A. S. Bradley, Phys. Rev. Lett. 119, 184502 (2017).
- (14) M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S. Bradley, Phys. Rev. Lett. 110, 104501 (2013).
- (15) D. R. Tilley and J. Tilley, Superfluidity and Superconductivity. (IoP, third edition, 1990).
- (16) C. F. Barenghi, L. Skrbek, and K. R Sreenivasan, Proc. Natl. Acad. Sci. U.S.A 111, 4647 (2014)
- (17) P. Walmsley, D. Zmeev, F. Pakpour, and A. Golov, Proc. Natl. Acad. Sci. U.S.A 111, 4691 (2014)
- (18) E. Fonda, K. R Sreenivasan, and D. P Lathrop, Proc. Natl. Acad. Sci. U.S.A 116, 1924 (2019)
- (19) A. W. Baggaley, C. F. Barenghi, and Y. A. Sergeev, Phys. Rev. E 89, 013002 (2014).
- (20) X. Rojas and J. P. Davis, Phys. Rev. B 91, 024503 (2015).
- (21) F. Souris, X. Rojas, P. H. Kim, and J. P. Davis, Phys. Rev. Appl. 7, 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. 124, 015301 (2020).
- (23) S. Forstner, Y. Sachkou, M. Woolley, G. I. Harris, X. He, W. P. Bowen, and C. G. Baker, New J. Phys. 21, 053029 (2019).
- (24) See appendix for details on experimental setup, pressure gradient and velocity calibration, discussion on two dimensionality, further analysis and extension of the theoretical model and additional data.
- (25) H. A. Kierstead, J. Low Temp. Phys. 23, 791 (1976).
- (26) D. J. Bishop and J. D. Reppy, Phys. Rev. Lett. 40, 1727 (1978).
- (27) S. Jose Benavides and A. Alexakis, J. Fluid Mech. 822, 364 (2017).
- (28) S. Babuin, E. Varga, L. Skrbek, E. Lévêque, and P.-E. Roche, EPL 106, 24006 (2014)
- (29) M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S. Bradley, Phys. Rev. Lett. 114, 155302 (2015).
- (30) A. Pouquet, R. Marino, P. D. Mininni, and D. Rosenberg, Phys. Fluids 29, 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. 168, 31-39 (2012).
- (32) K. W. Schwarz, Phys. Rev. B 31, 5782-5804 (1985).
- (33) R. J. Donnelly, Quantized Vortices in Helium II. Cambridge University Press, 1991.
- (34) R. J. Donnelly and C. F. Barenghi, J. Phys. Chem. Ref. Data., 27, 1217 (1998).
- (35) G. W. Stagg, N. G. Parker, and C. F. Barenghi, Phys. Rev. Lett. 118, 135301 (2017).
- (36) E. J. Hopfinger, F. K. Browand, and Y. Gagne, J. Fluid Mech. 125, 505 (1982).
- (37) J. Gao, W. Guo, S. Yui, M. Tsubota, and W. F. Vinen, Phys. Rev. B 97, 184518 (2018).
- (38) A. J. Groszek, M. J. Davis, D. M. Paganin, K. Helmerson, and T. P. Simula, Phys. Rev. Lett. 120, 034504 (2018).
- (39) R. N. Valani, A. J. Groszek, and T. P. Simula, New J. Phys. 20, 053038 (2018).
- (40) A. J. Groszek, T. P. Simula, D. M. Paganin, and K. Helmerson. Phys. Rev. A 93, 043614 (2016).
- (41) A. Cidrim, F. E.A. Dos Santos, L. Galantucci, V. S. Bagnato, and C. F. Barenghi, Phys. Rev. A 93, 033651 (2016).
- (42) E. Varga and L. Skrbek, Phys. Rev. B 97, 064507 (2018).
- (43) A. A. Berryman, Ecology 73, 1530 (1992).
- (44) A. J. Lotka, J. Phys. Chem. 14, 271 (1910).
- (45) H. Y. Shih, T. L. Hsieh, and N. Goldenfeld, Nat. Phys. 12, 245 (2016).
- (46) C. F. Barenghi, R. J. Donnelly, and W. F. Vinen, Phys. Fluids 28, 498 (1985)
- (47) M. L. Woyciekoski, L. A. M. Endres, A. V. de Paula and S. V.Möller Ocean Eng. 195, 106658 (2020)
- (48) R. M. B. Young and P. L. Read, Nat. Phys. 13, 1135 (2017).

## Supplementary Information

## Experimental setup and superfluid velocity calculation.

Flow in the Helmholtz resonators is driven and sensed capacitively using two aluminum electrodes deposited on the top and bottom wall of the device, forming a parallel plate capacitor (see Souris2017 for details on the fabrication process). An alternating voltage of amplitude applied to the electrodes of the device causes a periodic deformation of the walls of the basin due to electrostatic force, which pushes the superfluid in and out of the basin through the two side channels and into the bath, thus driving the Helmholtz mode. The Helmholtz resonance is observed as a periodic variation of the capacitance of the device.

The resonator is wired in a bridge circuit shown in Fig. SM1, balanced to the capacitance of the resonator at rest. Change in this capacitance caused by the Helmholtz resonance results in bridge imbalance and a current through the 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 . A standard 9V battery is used as the source of the bias voltage ( V). The capacitance bridge is the model General Radio 1615-A.

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 and the pressure gradient in the channels. Next, we calculate the superfluid velocity in the channels from the current through the detector.

Helmholtz equation of motion. – We derive the equation of motion of the superfluid in the Helmholtz resonator with explicit drive and damping forces. We approximate the flow in the channel as a mass on a spring, which is displaced by distance (positive in the direction away from the basin). The average displacement of the plates of the basin is denoted by (positive when the basin contracts). We begin by calculating the change in total density of the fluid inside the basin as a response to the mean deformation of the basin and the displacement of the superfluid inside the channel :

(5) |

Here is the mass of the fluid inside the basin, is the volume of the basin ( being its area and the confinement) and from the assumption that only the superfluid moves ( is the cross-sectional area of the channel; the factor of 2 comes from the two channels), and is the change in basin volume due to motion of the plates. A change in density corresponds to a change in pressure via the compressibility , , or

(6) |

Balancing forces on the plate (neglecting its inertia) yields

(7) |

where is the electrostatic force between the parallel plates of the capacitor formed by the circular electrodes in the basin, is the total applied voltage and N/m Souris2017 is the stiffness of the substrate deflection (note that this is double that in Ref. Souris2017 , where the stiffness refers to deflection of both plates in parallel). Expressing from Eq. 7 and substituting back in to Eq. 6 yields

(8) |

The superfluid inside the channel is accelerated by the pressure

(9) | ||||

(10) |

where we included a friction force . Here is a friction parameter with units of frequency but will remain otherwise unspecified for now. Rearranging,

(11) |

from which the resonance frequency follows

(12) |

where .

Finally, the driving pressure gradient, the quantity shown on the x-axis in Fig. 2A,B of the main text, is given by

(13) |

where we take only the component of the force on resonance with the Helmholtz mode ( being the AC drive), .

Calculation of velocity from detector current. – Whenever the capacitance bridge in the measurement circuit shown in Fig. SM1 becomes imbalanced, current will flow through the detector. Assuming that the bridge is tuned to the total capacitance of the devices at rest, the current through the detector at the frequency of the drive (the only component detected by the lock-in amplifier) is given only by the oscillation of the device capacitance due to the Helmholtz resonance,

(14) |

where is the bias voltage. The change of capacitance with superfluid displacement in the channel can be written as

(15) |

where is the capacitance of an undisturbed device with being the vacuum permittivity and dielectric constant of helium, respectively, and the area of the electrodes. The spacing between electrodes is given by , hence the second term in Eq. 15.

Neglecting the dependence of polarizability of helium on density Kierstead1976 and using the Clausius-Mossoti relation we can estimate the change in dielectric constant as

(16) |

The change of density with superfluid displacement can be calculated directly from Eq. 8 and using ,

(17) |

where . Differentiating Eq. 5 with respect to and putting the result equal to Eq. 17 results in an equation for and yields

(18) |

Inserting Eq. 18, Eq. 17 and Eq. 16 back into Eq. 15 yields

(19) |

Finally, the flow velocity is calculated as

(20) |

assuming that the background has been subtracted from .

### .1 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., is much larger than the coherence length, Å) for finite-size effects of 2D superfluidity to be relevant Bishop1978 . 2D turbulence, on the other hand, requires only that the fluctuating velocity is restricted to 2D. This is essentially controlled by the channel 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 Babuin2014 . Thus we first verify that the turbulence could be considered two-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 mm, which can, in principle, drive a 3D flow instability. It was shown by Benavides and Alexakis Benavides2017 , for systems of reduced dimensionality that the direction of the turbulent energy cascade critically depends on the ratio of forcing to confinement scale. Specifically, for the turbulence develops the 2D inverse energy cascade. Here, is the Reynolds number, which we define for our system using the effective quasi-classical viscosity He-II Babuin2014 as . In our experiments for nm and the highest experimentally achieved . Therefore, from a standpoint of classical turbulence, the turbulence in our devices ought to be in the 2D regime. It should be noted, however, that even if a few vertical modes of motion are possible, the inverse energy cascade responsible for appearance of large-scale features is still expected to be present Pouquet2017 .

The turbulent fluctuations, however, will be also strongly affected by the presence of quantized vortices, whose core size is on the scale of nm—significantly smaller than the confinement imposed by the device. A potential complication arises from vortex pinning on rough surfaces. The RMS surface roughness of our devices is expected to be Duh2012 about 1 nm, which puts the flow velocity required to dislodge a vortex from a typical surface defect Schwarz1985 at about 4 cm/s. The velocities we observe in the turbulent regime are significantly higher, thus it is unlikely that pinning plays an important role for our results.

It is in principle possible that a portion of the vortices in the flow are intrinsically three-dimensional, e.g. half-loops pinned on one of the opposing confining walls. To estimate the importance of such vortices we estimate their lifetime in a configuration shown in Fig. SM3(a). We assume a circular vortex attached to one wall aligned perpendicular to the applied oscillating flow. The self-induced velocity of the ring (neglecting pinning) as a function of its radius is given by Donnelly1991

(21) |

and, for stationary normal fluid, the change in radius is given by Donnelly1991

(22) |

where is the mutual friction constant Donnelly1998 and is the imposed superflow. We numerically integrate the evolution of for a range of initial radii and velocity amplitudes, terminating the calculation when either and the loop is annihilated or when m and the loop reconnects with the opposing wall, thus transforming into a vortex dipole. As shown in Fig. SM3(b), the typical lifetime of half-rings for the parameters typical of our experiment is shorter than the flow oscillation period , reaching, at most, about 0.6 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 feature in 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. This will cause the vortex to deform and be attracted to the opposing wall, thus slightly altering the lifetime. Changing the phase of the oscillating flow either does not significantly influence the outcome or causes loops of all sizes to quickly decay. Changing the angle between the plane of the loop and flow velocity would result in a somewhat more complicated transient flow, which is, however, unlikely to terminate in a significantly different manner.

The calculation outlined above is not valid for vortex loop radii comparable with the surface roughness nm, since the vortex will be subject to highly nonuniform flow resulting from the surface imperfections. Stagg et *al.* Stagg2017 studied a flow close to an irregular surface using a simulation of vortices in a Bose-Einstein condensate in the zero-temperature limit using the Gross-Pitaevskii equation (GPE). In that work, a dense layer of vortices was found in the rough landscape of the surface sustained by intrinsic nucleation of vortices on the protruding peaks of the surface. It is possible that such a dense boundary layer exists in our case as well, however, results obtained using GPE ought to be adopted with caution for helium at finite temperatures. Intrinsic nucleation of vortices in He-II requires significantly higher velocities and mutual friction at finite temperatures will strongly damp any small, highly-curved vortex structures. Regardless, this boundary layer is expected to be confined to within the scale of surface roughness Stagg2017 which in our case is about 0.1% of the confinement thus making it unlikely to be a significant contribution to the observed macroscopic drag.

Finally, the vortices connecting the two confining walls can, in principle, deform arbitrarily on the scale of the confinement, . We estimate the dynamical importance of these deformations by comparing their typical rate of decay to the time scale of their forcing, i.e., the flow oscillation period. Assume that the vertical modes of flow, mediated by the vortex deformation, take the form of a cascade of Kelvin waves—helical wave modes on vortices Donnelly1991 . The decay rate of a Kelvin wave mode of wave vector is Barenghi1985a . The smallest admissible results for K and nm in s. Increasing temperature will decrease . The decay rate of Kelvin waves is thus significantly faster than the time scale of their pumping (i.e., flow period, which is of the order of 1 ms) and comparable to the inverse frequency of the Kelvin mode itself Donnelly1991 , i.e., no Kelvin wave cascade is likely to develop along the individual vortices since the largest scales are already in the dissipative range. Other modes of vortex deformation (e.g., solitons Hopfinger1982 ), which cannot be decomposed to Kelvin waves, are possible. However, since the local velocity of the deformed line, and thus its decay rate mediated by mutual friction, are primarily determined by the local curvature, we expect the decay of these deformations to be comparable to that of the Kelvin waves. The amplitude of thermally excited Kelvin waves is also expected to be negligible Barenghi1985a . We therefore consider the vortices in our system which span the confinement to be point-like. Decreasing the temperature, particularly below 1 K, would suppress the mutual friction damping and allow the vortices to deform strongly. Therefore we expect the turbulence to cease to be 2D-like at sufficiently low temperatures.

### .2 Vortex density model

Discrete quantized vortices are transported by flow similar to how vorticity is transported in classical 2D flow Kraichnan1980 :

(23) |

where the terms on the right hand side correspond, respectively, to splitting, decay by collision and generation of vortices by the external drive. Note that this is not simply a passive scalar transport with source terms, since depends on the vortex distribution. We expect the splitting rate to depend on velocity, possibly exhibiting critical behavior itself. For simplicity, however, we take to be constant since in a high-velocity regime it is likely to be dominated by the flow oscillation period, which is independent of velocity.

Averaged over the flow oscillation period, the advection term will have no effect on the vortex density far from the system boundaries. In the region near the boundaries, however, some vortices will be transported toward the wall and annihilated, i.e., the average effect of the advection is to reduce the vortex number. Since the vortex density will vary on the scale of the channel width, we approximate the gradient term as . Putting , where characterises the inhomogeneous distribution of vortices throughout the channel, we recover Eqs. 1,2 of the main text. The assumption of velocity-independent and limits the applicability of the model to turbulent states at relatively high velocity and makes it unsuitable for modelling transition to turbulence from the laminar state or predicting the scaling of vortex number with velocity.

Following from Eqs. 1, 2 of the main text, the total vortex density and polarization obey

(24) |

and

(25) |

where we grouped all terms depending on to new terms and . The generation terms in Eqs. 1, 2 of the main text represent extrinsic or intrinsic nucleation of vortices and are likely to be concentrated near the sharp corners connecting the basin and the channel. The vortices generated at these edges in a polarized configuration are advected into the channel where they contribute to the observed drag. Near the corners, however, the polarization will likely be dominated by the instantaneous flow and, averaged over the flow oscillation period, making the , independent of . For simplicity we adopt and as independent control parameters, rather than . It should be noted, however, that it is the assumption of -independent that allows for bistable solutions.

The equations above are assumed to be local, but spatially averaged quantities are required for comparison with the experiment. The total vortex density is an always positive quantity and thus, to a first approximation, can be replaced by its spatial average. The vortex polarization , on the other hand, has a vanishing average since we assume that 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. The basin is removed and a single channel runs through the entire length of the device, but otherwise we assume general flow features similar to the real device (e.g., flow direction, behavior of the generation terms ). From the symmetry of the problem, is anti-symmetric with respect to mirroring about either of the axes and thus can be decomposed into orthogonal modes as

(26) |

and the spatial average is then given by . In the actual device geometry the modes in the expansion Eq. 26 will be more complicated but could, in principle, be constructed by a suitable transformation of the rectangular domain of Fig. SM4 onto the actual device geometry. Truncating the expansion at the lowest mode (which is likely to be the dominant term in the generation ) allows us to essentially use Eqs. 24, 25 as they are and recover the results from the main text.

In principle higher modes can be considered, where Eq. 25 would be replaced by a set of equations for each mode coupled through nonlinear terms. The generation term is unlikely to have a single-mode decomposition and the nonlinear terms (in ) in Eq. 25 will excite higher modes at the expense of lower modes. This picture is fully consistent with the forward enstrophy (quadratic integral of vorticity) cascade of classical 2D turbulence Kraichnan1980 . Higher modes will again exhibit near-degeneracy of the and solutions lifted by the appropriate mode of the and, possibly, by the lower-lying modes through the nonlinear terms. This will result in more general multi-stability of the mean vortex number , as illustrated in Fig. SM5b. As the flow velocity or drive increases, the system will randomly select either the or solution. As the drive increases further, the higher-order terms will become important, which will again be selected randomly, splitting each branch further. The beginning of this tree of turbulent states is, perhaps, already seen in the high-velocity part of the pressure-velocity curves at 1.4 K shown in Fig. 2a of the main text and highlighted in Fig. SM5 where three distinct turbulent branches are clearly seen.

### .3 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 with the 1067 nm confinement (shown in Fig. SM7 and Fig. 2 of the main text) under identical conditions. The bistability is again present in the 805 nm confinement, but in a weaker form and in the temperature range of 1.6–1.8 K. The bistability also tends to be suppressed at high velocities.

Within the model of Sec. .2 the bistability can be destroyed in several ways, apart from the already discussed temperature dependence of the decay parameter . For example, increasing the splitting rate above a critical value will result in a single solution with (i.e., frequent splitting will completely mix the flow). Similarly, increasing above a critical value will destabilise the solution with sign of opposite to the sign of (i.e., opposing polarization will be overwhelmed by the strong drive). Additionally, the confinement of the device likely affects the parameter as well – the smaller 805 nm device allows for lesser lateral deformation of the vortices and hence lowers the effective cross-section for collision, thus reducing , which tends to reduce the bistability.

In fact the bistability is not necessarily destroyed completely. If the environmental disturbances (e.g., vibration of the cryostat) are non-negligible compared to the relative stability of the less-stable state (controlled, for example, by the parameter), the flow would stochastically transition to the more stable state whenever a sufficiently strong fluctuation randomly occurs. This scenario is consistent with the fact that the transition between the two turbulent states in the temperature range 1.6–1.8 K for the 805 nm confinement does not appear to have a well defined critical velocity.

The device-dependence of the parameter discussed above, however, does not account for the complete lack of bistability at lower temperatures in the 805 nm device. One possibility for this observation is that critical velocity of type II (in Fig. 2(b) of the main text) moved beyond the critical velocity of type I. Once the laminar flow becomes unstable, only one turbulent state would be available which would thus be the only state observed. Indeed, this would be consistent with a relatively narrow hysteretic region at low temperatures in Fig. SM6. Additionally, Fig. 2(c) of the main text could suggest that the closing of the gap between critical velocities of types I and II is plausible even for the 1067 nm device at lower temperatures. However, due to lack of data from lower temperatures and lack of a model of the critical velocities this scenario ought to be regarded as a speculation at this point.

In order to describe the destruction of bistability precisely, significantly more detailed understanding of the critical velocities and the vortex-boundary interaction, and the parameters of the Eqs. 24, 25 that stem from it, would be required. This will depend, for example, strongly on the morphology of the surface Stagg2017 and is beyond the scope of this work.