Scaling regimes of active turbulence with external dissipation

Active fluids exhibit complex turbulent-like flows at low Reynolds number. Recent work predicted that 2d active nematic turbulence follows universal scaling laws. However, experimentally testing these predictions is conditioned by the coupling to the 3d environment. Here, we measure the spectrum of the kinetic energy, $E(q)$, in an active nematic film in contact with a passive oil layer. At small and intermediate scales, we find the scaling regimes $E(q)\sim q^{-4}$ and $E(q)\sim q^{-1}$, respectively, in agreement with the theoretical prediction for 2d active nematics. At large scales, however, we find a new scaling $E(q)\sim q$, which emerges when the dissipation is dominated by the 3d oil layer. In addition, we derive an explicit expression for the spectrum that spans all length scales, thus explaining and connecting the different scaling regimes. This allows us to fit the data and extract the length scale that controls the crossover to the new large-scale regime, which we tune by varying the oil viscosity. Overall, our work experimentally demonstrates the emergence of universal scaling laws in active turbulence, and it establishes how the spectrum is affected by external dissipation.

These issues have been the subject of intense research and debate 5,11,12,31,[34][35][36][37][38][39][40][41][42][43] . Inspired by bacterial turbulence, initial work on extended Toner-Tu equations for polar flocks showed that these fluids exhibit flow spectra with scaling regimes. However, the corresponding scaling laws have non-universal exponents, whose values depend on model parameters 11,36,37,[44][45][46][47] . In contrast, recent theoretical work * These authors contributed equally to this work † These authors contributed equally to the supervision of this work on active nematics predicted scaling laws with universal exponents, which are independent of the active fluid properties (e.g. viscosity or activity) 34,35 . Experimental evidence for such universal scaling, however, remained elusive, in part due to the challenge of accurately measuring flows for sufficiently broad ranges of scales.
In addition, most of the theoretical predictions and numerical studies have addressed the case of two-dimensional flows. As in classical turbulence, 2d systems are often embedded in a 3d setup. It is, therefore, necessary to separate the genuinely 2d features from the effects due to the coupling to the environment. In classical 2d turbulence, for instance, frictional dissipation with the environment not only cuts off the energy cascade towards large scales, but it also modifies the scaling exponent of the enstrophy cascade towards small scales 48,49 .
Here, we address the influence of external dissipation on the scaling regimes in microtubule-based active nematics powered by motor proteins 17,22 . In this setup, the 2d active nematic is surrounded by layers of oil and water, two passive fluids. Our comprehensive experimental study spans a broad range of spatial scales and includes systematic variation of the oil viscosity. We also present a theoretical framework that incorporates the hydrodynamic coupling with the environment and provides an explicit expression for the full spectrum of the turbulent active flows. Altogether our results provide experimental evidence of universal scaling laws in 2d active turbulence, and determine the ranges where these scaling laws are observable in terms of the physical parameters. Our predictions include six different scaling regimes, which we classify in terms of three length scales: the average vortex size, the height of a external fluid layer, and a viscous length that controls whether dissipation is dominated by either the active or the external fluid. This analysis reveals that, in contrast to classic turbulence, external dissipation does not just introduce a small-q cut off to the scaling behaviour, but it also yields a new scaling regime. In our experiments, we vary the oil viscosity over more than four orders of magnitude and observe three of the predicted scaling regimes. The remaining regimes might be observed either in alternative experimental setups or in other systems such as cell monolayers.
Beyond the scaling laws, the detailed comparison between theory and experimental data allows us to probe other open questions in active turbulence. Specifically, we seek to explain the emergence of a characteristic vortex size. The coupling to external fluids selects a characteristic wavelength at the onset of spontaneous flows 22,[50][51][52] . However, we show that this selection mechanism cannot fully account for the average vortex size observed in the turbulent regime. Thus, our analysis suggests that nonlinear effects in the active fluid, possibly involving energy transfer across scales, also contribute to vortex size selection. Moreover, we establish the range of validity of our theory, as we observe that it correctly fits the data for an intermediate range of oil viscosities but fails for extremely low and high viscosities. Our measurements suggest directions for future improvements to the theory, such as including vortex-vortex correlations.
In our experiments, we prepare an active nematic film by self-assembly of micrometer-long stabilized microtubules at the interface between a H water ≈ 40 µm-thick water layer on a glass slide and a H oil ≈ 3 mm-thick oil layer open to the air 53 (bottom inset in Fig. 1a, Methods). The microtubules are bundled under the depleting action of polyethylene-glycol (PEG), which facilitates crosslinking by kinesin molecular motors clustered with streptavidin (Methods). Fueled by adenosine triphosphate (ATP), the motors generate active shear stresses leading to extension and buckling of the microtubule bundles. The bundles then acquire local nematic order interrupted by half-integer topological defects, and the film exhibits disordered large-scale flows, which is the state of active nematic turbulence.
To study the statistical properties of active turbulent flows, we measure their so-called kinetic energy spectrum E(q) ∝ q |ṽ(q)| 2 , whereṽ(q) are the Fourier components of the flow field v(r), with wave vector q (Methods). To explore the different flow regimes, we vary an external parameter, namely, the oil viscosity, while keeping the intrinsic active fluid parameters -including motor, microtubule, and ATP concentrations -fixed. Therefore, the so-called active length, which compares the active and elastic stresses, is kept constant in all of our experiments. Upon increase of the oil viscosity, the entire kinetic energy spectrum decreases (Fig. 1), which is consistent with the previously-observed decrease of flow speed 53 . At low oil viscosities, the spectrum features at least three regimes: a large-scale (small-q) regime that is followed by a peak, an intermediate regime, and a crossover to a smallscale (large-q) regime (Fig. 1a). As oil viscosity increases, the peak shifts to smaller scales, expanding the range of the large-scale regime and shrinking the intermediate regime until it can no longer be observed for high oil viscosities (Figs. 1b and 1c). In parallel to these changes in the flow properties, we also observe a higher density of defects for higher oil viscosities (top insets in Fig. 1).
To understand our measured spectra, we develop a theoret-   Figure 2 Scaling regimes of turbulent flows in an active nematic film in contact with an external fluid. The different regimes are predicted at length scales (∼ 1/q) either larger or smaller than the mean vortex radius R * , the viscous length v = ηn/ηext, and the thickness H of the external fluid layer. This figure summarizes the scalings in the thick-layer limit qH 1; see Fig. S2 for the predictions in the thin-layer limit qH 1.
ical framework that accounts for the hydrodynamic coupling between the active film and the external water and oil layers, with their corresponding boundary conditions. Active flows in the nematic film induce flows in the passive layers which, in turn, influence the active-film flows (Fig. S1). Generalizing previous works 53-57 , we first obtain the Green function of the active flows subject to this feedback (SI Section A, see Eq. (S11)). We then establish a relationship between the velocity power spectrum in our coupled three-fluid system and the vorticity power spectrum of an isolated active nematic film (SI Section B, see Eq. (S22)). The latter was previously predicted by Giomi using a mean-field theory approach based on decomposing the vorticity field into a superposition of N uncorrelated vortices 34 . Based on simulation results, Giomi's theory assumes that each vortex has a uniform and size-independent vorticity ω v , and that vortex areas follow an exponential distribution with mean a * = πR 2 * , where R * is the mean vortex radius. For a configuration that has on average N vortices over a total system area A, we predict a kinetic energy spectrum (SI Section B) Here, I 0 and I 1 are modified Bessel functions, and B = N ω 2 v /(32π 3 A) is a prefactor related to the total enstrophy, which is independent of both the wave number q.
To discuss the scaling regimes predicted from Eq. (1), we consider a simplified situation with just one external fluid layer. The energy spectrum then depends on three length scales: the mean vortex radius R * , the layer thickness H, and the viscous length v = η n /η ext defined by the ratio of the two-dimensional viscosity of the nematic film, η n , and the three-dimensional viscosity of the external fluid, η ext . In the thick-layer limit, qH 1, we summarize the predicted scaling regimes in Fig. 2, which we discuss below. The thin-layer limit, qH 1, is discussed in Fig. S2 in SI Section C. Active flows are characterized by different scaling laws at scales smaller and larger than the mean vortex size R * (horizontal axis in Fig. 2). Furthermore, at length scales smaller than the viscous length v = η n /η ext , dissipation is dominated by the viscosity of the active film. As a result, in this regime, the scaling laws are those recently predicted and numerically demonstrated for isolated active nematic films 34,35 , with no effect of the external fluid (top half in Fig. 2). At length scales larger than v , however, dissipation in the external fluid dominates, yielding new scaling laws (bottom half in Fig. 2).
In our experiments, the water viscous length water = η n /η water is much larger than the largest length scale of our measurements: q water 1 for all q values. Therefore, the flows in the water layer do not produce any new scaling regimes in our experiments. In contrast, varying oil viscosity over orders of magnitude, we reach oil viscous lengths oil = η n /η oil that fall within our measurement window. Consequently, the flows in the oil layer are responsible for some of the scaling regimes that we observe. Our measurements probe length scales that are smaller than the oil layer thickness H oil ≈ 3 mm. Thus, our experiments operate in the thick-layer limit, with scaling properties as predicted in Fig. 2.
Consistent with these predictions, we experimentally observe E(q) ∼ q −4 at small scales, and E(q) ∼ q 1 at large scales, for all oil viscosities (Figs. 3a and 3b). For low and intermediate oil viscosities, we also observe signatures of the E(q) ∼ q −1 regime at intermediate length scales, larger than the vortex size R * but smaller than the viscous length oil (Fig. 3c). The q −4 and q −1 scaling behaviors are intrinsic properties of an active nematic film, respectively characterizing the small-scale flows inside vortices and the large-scale flows due to hydrodynamic interactions in the film 35 . In contrast, the q 1 scaling stems from the hydrodynamic coupling to an external fluid. All these scaling relations involve universal exponents, which are independent of the properties of the fluids. Consistently, by varying the oil viscosity, we tune the range of the different regimes without changing their scaling exponents.
Having demonstrated the scaling regimes, we quantitatively compare our prediction for the full energy spectrum (Eq. (1)) to the experimental data. Knowing the values of η water , η oil , H water , and H oil , we use B, η n , and R * in Eq. (1) as fitting parameters. Despite its assumptions, our theory fits the data remarkably well for a wide range of intermediate oil viscosities (9.7 × 10 −3 Pa·s < η oil < 0.39 Pa·s, see Fig. S3), as exemplified in Fig. 3c. Visual inspection of the experiments suggests smaller vortices as η oil increases (see insets in Fig. 1 and Movie S1). Yet, the mean vortex radius obtained from these fits, R s * , is independent of oil viscosity (Fig. 3d). Finally, the oil viscous length oil = η n /η oil decreases slightly faster than ∼ 1/η oil (Fig. 3e), indicating that η n weakly decreases with oil viscosity (Fig. S4). The physical origin of this dependence remains an open question.
To try to rationalize the fact that R s * is independent of oil viscosity, we perform a linear stability analysis. Due to the coupling to the external fluid, the linear growth rate of the spontaneous-flow instability that powers active nematic turbulence acquires a maximum at finite wavelengths 50-52 (SI Section D and Fig. S5a). This maximum selects a characteristic scale λ m of flow patterns at the onset of turbulence 22 (Fig. S5b). This linear analysis predicts that λ m changes with oil viscosity, in contrast with experimental observations for R * . However, nonlinear effects may modify this selected scale at later stages. In stationary fully-developed turbulence, earlier work 35 has shown that vortex size is given by the critical wavelength of the instability, λ c . We now show that this length scale is largely independent of oil viscosity (SI Section D and Fig. S5c).
To address this issue experimentally, we directly measure the distributions of vortex areas n(a) 19,21,27 (Fig. 4a, Methods). Fitting their exponential tails as n(a) ∝ exp(−a/a * ) 34 , we obtain a mean vortex radius R v * = a * /π as a function of oil viscosity (Fig. 4b). The results show a decrease of vor-tex radius for high oil viscosities (Fig. 4b), outside the range of validity of the presented model. For lower oil viscosities, within the range of validity of the model, the trend is weaker, suggesting that nonlinear effects contribute to vortex size selection in our system. Our analysis thus calls for further research to understand vortex size selection in active turbulence.
Finally, we use our measurements to examine some of the assumptions of the theoretical framework. We observe that vortex area distributions have exponential tails (Fig. 4a), in agreement with both the theoretical assumption and previous experiments 19,21,27 . Here, we find that this feature does not change with the oil viscosity. We also measure the correlation functions for velocity and vorticity 39,58,59 (Figs. 4c and 4d) and obtain the corresponding correlation lengths (insets in Figs. 4c and 4d), which exhibit dependencies on the oil viscosity that are very similar to that of the vortex size R v * (Fig. 4b). This observation suggests that these lengths are all proportional to one another (Fig. S7), validating our theoretical assumption that there is a unique scale related to activity (R * ). However, at intermediate and high oil viscosities, we measure negative correlations of the vorticity field at distances comparable to, and even larger than, the vortex size (Fig. 4d). Thus, our data calls for future work on the theory of active turbulence, going beyond mean-field approximations and accounting for vortexvortex correlations. Other improvements to the theory could take into account the structure of defect cores, and include a more accurate treatment of flow alignment effects.
In summary, we have experimentally measured and theoretically explained universal scaling laws in active turbulence, thus drawing parallels to classical turbulence. Specifically, we have found scaling regimes that are intrinsic to an active nematic film as well as other regimes that result from the coupling to an external fluid. In addition, we have developed a theoretical framework that provides the rationale for the different regimes and yields an explicit form for the turbulent spectra as a function of parameters. By fitting the predictions to the data, we have extracted the crossover length scales, and shown how they depend on the viscosity of the external fluid. Our analysis of these results paves the way toward addressing open questions in active turbulence, from vortex size selection to the role of vortex-vortex correlations, and to pursue a deeper understanding of the fundamental similarities and differences between inertial and active turbulence.

COMPETING INTERESTS
The authors declare no competing interests.

DATA AVAILABILITY
All data are available from the authors upon request.

CODE AVAILABILITY
All codes are available from the authors upon request.

Preparation of the active nematic
Clusters of biotinylated kinesin-1 (K401-BCCP-6His) were assembled with tetrametric streptavidin at a ∼2:1 ratio. Afterwards, these molecular motors were mixed with adenosine triphosphate (ATP) and an ATP-regenerating system (Pyruvate Kinase/Lactic Dehydrogenase enzymes (PK/LDH) and phosphoenol pyruvate (PEP)). The non-adsorbing polymer polyethylene glycol (PEG) (20 kDa) was used as depleting agent. To avoid photobleaching and protein oxidation, a mixture of antioxidants (glucose, catalase, 1,4-dithiothreitol (DTT), Trolox and glucose oxidase) was also included. Biocompatibility of the active nematic layer with the oil interface was assured with the surfactant Pluronic F-127. This final suspension was mixed with ∼1.5 µm microtubules stabilized with guanosine-5'-[(α,β)-methyleno]triphosphate (GM-PCPP) (Jena Bioscience, NU-405S), of which the 0.8% was labeled with the Alexa-647 fluorophore. The scattered fluorescent microtubules formed a speckle pattern, from which we measure the velocity field through particle image velocimetry (PIV). Final compound concentrations are listed in Table I.
The active nematic (AN) was finally assembled by depositing a 1.5 µL droplet onto a polyacrylamide-functionalized glass slide within a 10 mm wide and 10 mm high polypropylene cylinder, and covered with 300 µL of polydimethylsiloxane oil (PDMS) with a viscosity in the range of 6.3·10 -4 Pa·s to 12 Pa·s. Micotubules spontaneously adsorbed onto the oil/water interface leading to the formation of the AN. The mean height of the water layer was obtained considering a cylinder and measuring the cross section area of the drop through fluorescence microscopy images. In the case of the oil layer, its height was directly measured with a millimetric ruler. We took the height at the center, where the AN drop is placed, neglecting the meniscus.
To have intermediate oil viscosities, oil mixtures were prepared. The final viscosity was estimated using the Arrhenius mixing rule log(η 12 ) = x 1 log(η 1 ) + (1 − x 1 ) log(η 2 ), where η 12 , η 1 and η 2 are the oil viscosities of the oil mixture and of the mixed compounds 1 and 2, respectively, and x 1 is the molar fraction of oil 1 60,61 .

Observation of the active nematic
The active nematic layer was imaged by means of fluorescence microscopy (Nikon Eclipse Ti2-U) with an Andor Zyla 4.2 Plus camera controlled with the open-source software ImageJ Micro-Manager 62 . As light source, we used a red LED coupled to a Cy5 cube filter. Images were acquired typically at a frame rate of 2 fps and with a spatial resolution of 2.14 µm/px. In the case of experiments with high oil viscosities (>1 Pa·s) the frame rate was decreased to 1 fps and the spatial resolution was increased to 1.28 µm/px.

Data analysis
Raw images were treated with the open source software ImageJ 63 . In general, light intensity was equalized by dividing the intensity of each frame by the time average of the sequence. Further noise removal was achieved with a mean filter with a width of 2-4 px. Afterwards, the velocity field was computed with PIVlab for Matlab 64 . PIV window size was set as 1/64 th of the lateral system size. To reduce noise, we used the option 5 x repeated correlation. The velocity field was finally smoothed using the smoothn function developed by Garcia 65 . The angle-averaged spectrum of the kinetic energy density per unit mass, E(q), was computed from the obtained velocity fields as: with |ṽ(q)| 2 = 1 2π ϕ |ṽ(q, ϕ)| 2 ∆ϕ, where L 2 is the area of the field of view, · indicates a time average,ṽ are the Fourier components of the velocity v, and q and ϕ are, respectively, the magnitude and azimuth of the wave vector q. Note that the experimental Fourier transform is a discrete Fourier transform, whereas we use a continuous Fourier transform in the theory. This explains the difference in the prefactors between Eq. (2) and Eq. (S15) in the Supplementary Information. For each experiment, we averaged a total of 500 frames. Afterwards, Eq. (1) was fitted to the experimental E(q) with Mathematica v10.
Identification and characterization of vortices were carried out using the Okubo-Weiss (OW) parameter, as described previously 19,21,34 . Briefly, OW was calculated as Regions with OW < 0 were vortex candidates. To determine whether such regions were vortices, we checked if they featured a singularity by computing the winding number W (r) = 1/2π 2π 0 atan(v y (r, ϕ)/v x (r, ϕ))dϕ of the velocity field at a distance r = 24 px from their center. If a region with OW < 0 had a winding number W ∈ [0.95, 1.05], it was accepted as a vortex. The total area of each swirl was determined by the connected area with OW < 0. This vortex locator algorithm, allowed us to extract a distribution of vortex areas n(a) = N (a)/ a N (a), where N (a) is the number of vortices with area a. Since n(a) follows an exponential distribution, we extracted a characteristic vortex area a * and radius R v * = a * /π.
Supporting Information for "Scaling regimes of active turbulence with external dissipation" SUPPLEMENTARY TEXT

A. Hydrodynamic Green's function
In our experimental setup, a thin film of the active nematic fluid is in contact with a thicker layer of water underneath, and a much thicker layer of oil above. The water layer is supported by a solid substrate, and the oil layer is in contact with air (Fig. S1). Active flows in the nematic film induce flows in both passive fluid layers, which in turn, influence the flows in the active film. To account for this hydrodynamic coupling, in this section we obtain the Green function for the flow field in a two-dimensional viscous fluid film in contact with three-dimensional layers of other viscous fluids. This calculation generalizes previous work 57 by considering the simultaneous hydrodynamic coupling to two fluid layers with different boundary conditions, as in our experimental setup (Fig. S1).
All flows in our system take place at very low Reynolds numbers. Therefore, inertial forces are negligible (Stokes limit), and momentum conservation reduces to force balance. For the active nematic film, with two-dimensional shear viscosity η n and incompressible flow field v(r), force balance can be written as (S1) Here, the first term accounts for the viscous stress within the nematic film, P is the film's two-dimensional pressure, and f water and f oil are the viscous force densities exerted by the water and oil layers, respectively, on the active fluid film. Finally, f is the remaining force surface density acting on the film. In our system, this force density results from stresses associated with the orientational order of the nematic film, including elastic, flow-alignment, and active stresses. Here, we derive the Green function treating f as a source, without specifying its expression. Rather, we express the flow field of the nematic film in terms of the source force f as where G αβ (r−r ) is a hydrodynamic Green function, namely a generalization of Oseen's tensor.
Greek indices indicate spatial components, and summation over repeated indices is implicit. In Fourier space, we havẽ where we have introduced the Fourier decomposition as To obtain the Green function, we must obtain the viscous stress exerted by the water and oil layers on the active fluid respectively. Here, u(r, z) is the three-dimensional flow field of the passive fluids, and the subscript indicates the components along the active film's plane (z = 0, see Fig. S1). The viscous flows in the passive layers obey the Stokes equation with p the three-dimensional pressure. These flows are driven by the hydrodynamic coupling with the active film: u (r, 0 + ) = u (r, 0 − ) = v(r). As shown in Refs. 55,56 , if the film's flow is incompressible, ∇ · v = 0, then the pressure in the fluid layers is uniform, ∇p = 0, and the out-of-plane component of the velocity vanishes everywhere, u z = 0. Therefore, the layers' flow field is planar and harmonic; it obeys Laplace's equation: To obtain the Green function in Fourier space, as in Eq. (S3), we solve Eq. (S7) in terms of the planar Fourier modes of the flow field,ũ(q, z), which obey The water layer is in contact with a solid substrate at z = −H water , where we assume a no-slip boundary condition: u (r, −H water ) = 0. Respectively, the oil layer is in contact with air at z = H oil , where we assume a no-shearstress boundary condition: ∂ z u (r, z)| z=Hoil = 0. With these boundary conditions and u (r, 0 + ) = u (r, 0 − ) = v(r), the solutions to Eq. (S8) arẽ These solutions show that the flow penetrates into the oil and water layers to depths given by the inverse of the in-plane wave number, 1/q, unless limited by the layers' thicknesses, H water and H oil , as illustrated in Fig. S1. Introducing Eq. (S9) into Eq. (S5), we obtaiñ f water (q) = −η water q coth(qH water )ṽ(q), (S10a) f oil (q) = −η oil q tanh(qH oil )ṽ(q). (S10b) Finally, introducing these results into the Fourier transform of Eq. (S1), and using the incompressibility condition q ·ṽ = 0, we obtain the Green functioñ G αβ (q) = δ αβ − q α q β /q 2 η n q 2 + η oil q tanh(qH oil ) + η water q coth(qH water ) .
(S11) This function describes the hydrodynamic interactions in the active nematic film, accounting both for the direct interactions due to its incompressibility and viscosity, and also for the indirect interactions mediated by the the oil and water layers. The ratios between the viscosity of these external fluids and the two-dimensional viscosity of the nematic film define two viscous length scales: oil = η n /η oil , water = η n /η water . (S12) At scales larger than this viscous length, dissipation in the external fluid layer (either oil or water) dominates over dissipation within the nematic film.

B. Kinetic energy spectrum
In this section, we derive an analytical expression for the kinetic energy spectrum of the turbulent flows in our system. To this end, we establish a relationship between the flow field of our coupled three-fluid system (Fig. S1) and the vorticity field of an isolated active nematic film. By means of this relationship, we combine the hydrodynamic Green function obtained in SI Section A with Giomi's mean-field theory of active nematic turbulence 34 to predict the flow power spectrum in our experimental system. The kinetic energy per unit mass density, E, of the active nematic flows is given by Using the Fourier modesṽ(q) of the flow field, as introduced in Eq. (S4), the angle-averaged spectrum E(q), with q = |q| the wave number, is defined by where A is the area of the system, and · averages over realizations. In states where correlations of the flow field are isotropic, E(q) is given by To obtain the velocity power spectrum, we use Eqs. (S3) and (S11), which give Here, we have introduced the notation Λ(q) ≡ η n q 2 + η oil q tanh(qH oil ) + η water q coth(qH water ).
(S17) To eliminate the source force density f , and thereby obtain a closed-form expression for the velocity power spectrum, we leverage the force-balance condition for the active nematic film alone, without external fluids: (S18) Here, the subscript i indicates that the active nematic film is isolated. As in Eq. (S1), f accounts for source force density due to elastic, flow-alignment, and active stresses, albeit now in the isolated film. The source force of the isolated film (in Eq. (S18)) coincides with that of the full problem with external fluid layers (in Eq. (S1)) only if we ignore flow-alignment effects. In this limit, f does not depend explicitly on the flow field, and hence it is insensitive to the presence of external fluid layers. To exploit this fact and therefore be able to derive a closed form for the velocity power spectrum, we first ignore the flow alignment coupling. Following Ref. 35 , we take the curl of Eq. (S18) to obtain a Poisson equation for the vorticity field ω =ẑ · (∇ × v): where s is the vorticity source due to nematic forces. In Fourier space, this equation is written as Hence, the vorticity spectrum of the isolated active film is given by Comparing to Eq. (S16), we obtain with Λ(q) given by Eq. (S17). As explained above, this result is exact in the limit of vanishing flow alignment coupling.
In the presence of flow alignment, Eq. (S22) defines a closed approximation whereby the nematic forces, both passive and active, are included exactly, while the flow alignment forces are approximated by the hydrodynamics of the isolated problem. Equation (S22) relates the velocity spectrum of the active nematic film coupled to the external fluid layers with the vorticity spectrum of an isolated active nematic film. For the latter, we can now make use of the mean-field theory introduced by Giomi 34 , which is based on decomposing the vorticity field into a superposition of N uncorrelated vortices. Based on simulation results, the theory assumes that each vortex has a vorticity ω v independent of its size, and that vortex areas follow an exponential distribution with mean a * = πR 2 * , where R * is the mean vortex radius. With these assumptions, Giomi's mean-field theory predicts 34 where I 0 and I 1 are modified Bessel functions of the first kind. Introducing this result into Eq. (S22), and into Eq. (S15), we obtain the kinetic energy spectrum ) is a prefactor related to the total enstrophy, and independent of both the wave number q and the mean vortex radius R * .   Figure S2 Scaling regimes of turbulent flows in an active nematic film in contact with an external fluid. The different regimes are predicted at length scales (∼ 1/q) either larger or smaller than the mean vortex radius R * , the viscous length v = ηn/ηext, and the thickness H of the external fluid layer. This figure summarizes the scalings in the thin-layer limit qH 1; see the main text and Fig. 2 for the predictions in the thick-layer limit qH 1.

C. Scaling laws
In this section, we extend the discussion of the predicted scaling regimes given in the main text. As in the main text, we consider a situation with just one external fluid layer and classify the scaling regimes in terms of three characteristic lengths: the mean vortex radius R * , the viscous length v = η n /η ext , and the external layer thickness H. In the main text, we discuss the scaling regimes at scales much smaller than the layer thickness, qH 1, summarizing our results in Fig. 2. Here, we give the results in the opposite limit, looking at scales much larger than the layer thickness, qH 1. In this thin-layer limit, the scaling laws depend on the boundary condition of the external fluid at the non-active interface. We considered both a no-slip boundary condition, as for the water-substrate interface in our experiments, and a free surface, as for the oil-air interface in our experiments (Fig. S1). We summarize the results for all these situations in Fig. S2, which shows two additional scaling regimes, E(q) ∼ q 3 and E(q) ∼ q 0 , for no-slip boundary conditions.

D. Vortex size selection
In this section, we discuss the physical origin of the mean vortex size R * , a key parameter in the description of active turbulent flows. Active nematic turbulence results from the well-known spontaneous-flow instability in active nematic fluids 1,2,4 . Here, we ask whether vortex size is selected at the early stages of the instability, i.e. by its linear dynamics, or later on by the nonlinear dynamics.
In isolated active nematic films, the spontaneous-flow instability is a long-wavelength instability, whose growth rate is maximal at the longest wavelengths. As a result, the linear regime of the instability does not select any intrinsic wavelength but rather produces flow patterns with wavelengths determined by the system size. Indeed, system-sizedependent, stationary vortex patterns were observed in simulations at small activities, just past the instability threshold 35 . At high activity, however, these stationary vortices were unstable and evolved into turbulent flows. Simulations without flow-alignment revealed a sequence of instabilities whereby vortices break into smaller vortices, down to a characteristic vortex size proportional to the critical wavelength of the instability, λ c , which is itself proportional to the active length a 35 . In the presence of flow alignment, a similar transient cascade may exist. Moreover, flow alignment introduces a new nonlinearity in the problem enabling an additional potential mechanism of energy transfer in the steady state. It is thus reasonable to expect that, also in the presence of flow alignment, the vortex size for an isolated nematic is selected by a nonlinear mechanism.
This scenario is modified when the active nematic is not isolated but coupled either to a frictional substrate 58,66 or to external fluids 22,[50][51][52] , as in our system. The external fluids are isotropic and in the Stokes regime, and therefore they cannot transfer energy across scales at the steady state. However, the coupling to an external fluid modifies the spontaneous-flow instability of the active nematic at the linear level, making the growth rate achieve a maximum at finite wavelengths 51,52 . Consequently, an intrinsic wavelength λ m is selected from the very onset of spontaneous flows, which dominates the early stages of turbulence development 22 . To what extent nonlinear effects may modify this characteristic scale at later stages of turbulence is an open question that we can now analyze in light of our data.
To this end, we derive the wavelength λ m that is selected in the linear regime by the coupling between the active nematic and the external fluids in our experiments. We then compare these predictions to our experimental measurements of the vortex size in fully-developed turbulence at different oil viscosities.
To obtain the selected wavelength, we analyze the dynamics around the uniformly-oriented quiescent state of the active nematic film. Ignoring flow alignment and topological defects, the dynamics of the angle θ of the nematic director field n = (cos θ, sin θ) can be written as 35 The right-hand side of this equation accounts for the relaxation of distortions of the director field, which generate elastic nematic stress. We describe this stress in the approximation of one Frank constant K which, together with the rotational viscosity γ, controls the director relaxation rate 67 . Respectively, the second and third terms on the left-hand side of Eq. (S25) account for the advection and corotation of the nematic director by the flow field v, with vorticity ω =ẑ · (∇ × v).
In turn, the flow field is driven by nematic forces f as spec-ified by the force balance Eq. (S1). Ignoring flow alignment as in Eq. (S25), the nematic forces are given by The first contribution accounts for elastic nematic stresses described by the antisymmetric part of the the stress tensor 67 , Here, h = −δF n /δn = K∇ 2n is the molecular field computed from the Frank free energy for nematic elasticity which, in the one-constant approximation, is given by 67 Respectively, the second contribution in Eq. (S26) corresponds to the active nematic stress 4,68-70 where ζ is the active stress coefficient, and q αβ = n α n β − 1/2 δ αβ is the nematic orientation tensor in two dimensions, with cartesian components To perform a linear stability analysis, we introduce perturbations around the reference state as θ = 0+δθ and v = 0+δv, and we decompose them into Fourier-Laplace modes as δθ(r, t) = dΩ 2π δv(r, t) = dΩ 2π with wave vector q and frequency Ω. In terms of these modes, and to first order in perturbations, Eq. (S25) is written as In turn, the Fourier modes of the flow field are related to those of the nematic force,f , via Eqs. (S3) and (S11). To first order in perturbations, Eqs. (S26) to (S30) yield Introducing these results into Eq. (S3), using Eq. (S11), and then introducing the resulting δṽ into Eq. (S32), we obtain the growth rate of perturbations: where φ is the angle formed by the wave vector q and the directorn, such that q ·n = q cos φ, and Λ(q) is the hydrodynamic kernel given by Eq. (S17).
To get insight into wavelength selection, we consider a situation with only one external fluid layer with thickness H → ∞ and viscous length v = η n /η ext . In this situation, Eq. (S34) takes the simpler form Rescaling length and time by the active length and time, respectively, the growth rate Eq. (S35) can be expressed in dimensionless form as Here,Ω = Ωτ a , andq = q a are dimensionless variables. We have also introduced three dimensionless parameters: the viscosity ratio r ≡ γ/η n , the dimensionless viscous length v = v / a , and the sign of the active stress sgn(ζ) = ±1 for extensile and contractile stresses, respectively.
The direction of most unstable perturbations is along the director (φ * = 0) for extensile stresses (ζ > 0), and perpendicular to the director (φ * = π/2) for contractile stresses (ζ < 0). Along the most unstable direction, the growth rate Eq. (S37) has the shape shown in Fig. S5a. This figure shows that the coupling of the active nematic to the external fluid layer, represented by a finite¯ v = v / a , produces a maximum of the growth rate at a finite wavelength. This selected wavelength, λ m = 2π/q m , is obtained from the single real solution of the following cubic equation: As shown in Fig. S5b, and also apparent in Fig. S5a, the selected wavelength λ m has a non-monotonic dependence on the viscous length v . In the common situation in which the viscous length is larger than the active length, the selected wavelength decreases with the viscosity of the external fluid. However, when the viscous length becomes smaller than the active length, the selected wavelength increases with the viscosity of the external fluid. The coupling to the external fluid not only gives rise to the selected wavelength λ m but also modifies the critical wavelength of the instability, λ c = 2π/q c . This wavelength is determined by the condition Ω(q c ) = 0 along the most unstable direction, which gives with a given in Eq. (S36). As shown in Fig. S5c, the critical wavelength λ c increases monotonically with the viscosity η ext of the external fluid. At small η ext , i.e. large v = η n /η ext , the critical wavelength saturates at its value for isolated nematic films: lim v a λ c = 2π 2 + r/2 a . (S40) In the opposite limit, when dissipation is dominated by the viscosity of the external fluid, the critical wavelength becomes proportional to the external fluid viscosity: In summary, the coupling to the external fluid endows the spontaneous flow instability with a linear wavelengthselection mechanism. This mechanism selects the characteristic scale of turbulent flows at early stages 22 . However, nonlinear effects may modify the selected scale at later stages, and thereby dictate the vortex size in the stationary, fullydeveloped turbulent regime. To assess this point in our experiments, we analyze how stationary vortex size varies with oil viscosity.
If stationary vortex size were dictated by the linear wavelength-selection mechanism, vortex size should vary with oil viscosity as in Fig. S5b. In our experiments, the oil viscous length is oil R * , and the vortex size should be R * a . Therefore, we are in the regime v > a , and vortex size should moderately decrease with oil viscosity (Fig. S5b). Instead, if vortex size were dictated by the nonlinear selection mechanism of Ref. 35 , vortex size should be proportional to the critical wavelength, and therefore vary with oil viscosity as in Fig. S5c. Thus, again in the regime v > a , vortex size should remain roughly independent of oil viscosity (Fig. S5c).
Our experimental measurements show that the mean vortex radius, as measured from the vortex area distribution, is rather independent of oil viscosity over a range spanning several orders of magnitude of oil viscosity (Fig. 4b). At very high oil viscosities, we observe a moderate decrease of vortex radius (Fig. 4b), whose interpretation remains an open question. This analysis suggests that the linear selection mechanism alone does not explain the observed vortex size in the stationary regime, thus implying some form of energy transfer across scales, be it transient and/or steady. Elucidating the mechanism of vortex size selection in more detail thus remits to fundamental open questions in active turbulence that we defer to future work.  Figure S5 Vortex size selection. a, Growth rate (Eq. (S37)) along the most unstable direction for the spontaneous-flow instability in an active nematic film coupled to an external fluid layer. The viscosity ratio is set to r = γ/ηn = 1, and flow alignment is ignored. The different curves correspond to different values of the external fluid viscosity ηext, expressed in terms of the viscous length v = ηn/ηext. The blue line corresponds to an isolated active nematic film, without external fluid. Lengths and time are rescaled by the active length and time, a and τa defined in Eq. (S36). b, Wavelength at which the growth rate is maximum, as obtained from Eq. (S38), as a function of the viscous length. This wavelength is selected by the linear dynamics upon the spontaneous-flow instability. c, Critical wavelength of the spontaneous-flow instability (Eq. (S39)) as a function of the viscous length. Figure S6 Comparing vortex sizes. Ratio of the mean vortex radii obtained from the spectral fits (R s * , Fig. 3d) and from the vortex area distribution fits (R v * , Fig. 4b). In the range of oil viscosities in which the theory fits the data well, both estimates of R * exhibit similar trends with oil viscosity. Hence, their ratio is rather independent of oil viscosity. Figure S7 Comparison of characteristic length scales of active nematic turbulence. The different length scales are defined as follows: R s * is the mean vortex radius obtained by fitting the kinetic energy spectra (Figs. 3c and 3d); R v * is the mean vortex radius obtained by fitting the exponential tail of the vortex area distribution (Figs. 4a and 4b); vv is the distance at which the velocity autocorrelation is 0.5 (Fig. 4c); ωω is the distance at which the vorticity autocorrelation is 0.5 (Fig. 4d); q −1 max is the inverse of the wave number for which the kinetic energy spectrum has its maximum; oil = ηn/ηoil is the oil viscous length obtained from the fits of the kinetic energy spectra (Fig. 3c). Except for oil, all other length scales have similar behaviours with oil viscosity, consistent with them being proportional to one another.
Movie S1 Flow field of the AN at different oil viscosities. As the oil viscosity is increased, the flows within the active layer become slower and the characteristic vortex size decreases. The velocity fields are the ones obtained from the PIV analysis of the experiments. Colors of the vectors indicate their magnitude. For the sake of a better visualization, we show a field of view smaller than the one used to compute the kinetic energy spectra and plot one vector every three. The movie is sped up x5.