Active hydraulics laws from frustration principles

Viscous flows are laminar and deterministic. Robust linear laws accurately predict their streamlines in structures as complex as blood vessels, porous media and pipe networks. However, biological and synthetic active fluids defy these fundamental laws. Irrespective of their microscopic origin, confined active flows are intrinsically bistable, and therefore non-linear. As a consequence, their emergent patterns in channel networks are out of reach of available theories, and lack quantitative experiments. Here, we lay out the basic laws of active hydraulics. We show that active hydraulic flows are non-deterministic and yield degenerate streamline patterns ruled by frustration at nodes with an odd coordination number. More precisely, colloidal-roller experiments in trivalent networks reveal how active-hydraulic flows realize dynamical spin ices. The resulting streamline patterns split into two distinct classes of self-similar loops, which reflect the fractionalization of topological defects at the subchannel scales. Informed by our measurements, we formulate the laws of active hydraulics as a double spin model. A series of mappings on loop O(n) models then allow us to exactly predict the geometry of the degenerate streamlines. We expect our fundamental understanding to provide robust design rules for active microfluidic devices, and to offer unanticipated avenues to understand the motion of living cells and organisms in complex habitats.

When they invented irrigation, the ancient civilizations of Egypt and Mesopotamia relied on the first fundamental rule of hydraulics: mass conservation [1].In our modern language we formulate it as the first Kirchhoff's law.In steady state, the sum of the fluxes vanishes at each node of a channel network ( j Φ ij = 0, where Φ ij is the flux from node j to node i).It is only eight millennia later, that Hagen, Poiseuille and Darcy discovered the second law of hydraulics [2][3][4].Darcy's law relates the mass fluxes to pressure gradients: Φ ij = −K ij (P i − P j ), where P i is the fluid pressure and K ij the hydraulic conductance.Given these two linear laws, and a set of boundary conditions, we can predict the viscous flows of any hydraulic network, regardless of its geometrical complexity.Although it can be computationally challenging to solve and optimize, viscous hydraulics is deterministic and fully predictable.
However, over the past decade physicists, chemists and material scientists have complexified this neat picture.They engineered active fluids that escapes the fundamental laws of hydraulics [5].Either by hacking biological engines, such as molecular motors and bacteria, or by motorizing synthetic soft materials, we have learned how to power fluids at the scale of their elementary building blocks.When confined, the resulting active materials enjoy spontaneous laminar flows even in the absence of any external drive or boundary motion [3,[7][8][9][10][11][12].In channels and pipes, active fluids flow in one direction or the other with the same probability, and can even resist opposing pressure gradients [13].In other words, they are intrinsically bistable [9,13,14].This fundamental deviation from Darcy's law makes active hydraulics a challenging problem whose basic laws are yet to be determined.
From an experimental perspective, active hydraulics has been mostly restricted to straight channel geometries and closed loops [7,9,15].To the best of our knowledge, the investigation of interconnected channel networks has remained limited to pioneering experimental demonstrations in simple one-node networks [9].More significant progress have been made on the theory front, in particular by Dunkel, Woodhouse and coworkers, who realized the computational power of Active Fluidic Networks (AFN) and their dynamical multistability [14,[16][17][18].
In this article, we perform large-scale active hydraulics experiments.We show that spontaneous laminar flows are frustrated in networks including nodes with an odd coordination number.Focusing on fully frustrated networks, we show that the resulting active flows realize dynamical spin ices with extensively degenerate flow patterns.We show that unlike passive fluids, the random streamline geometry depends on the aspect ratio of the elementary channels.We explain this polymorphism by combining experiments and numerical simulations, and show that it originates from topological-defect fractionalization at the subchannel level.We then elucidate the self-similar geometry of the flow patterns by mapping them on the frustrated structures of magnetic spin ices, and of the so-called loop O(n) models [8,9,20,21].Altogether our findings allow us to identify the full set of laws ruling the steady flows of active matter circulating through interconnected channels networks.
To tackle the problem of active hydraulics experimentally, we need a model active fluid, and model hydraulic networks, see Fig. 1.We provide all the technical details about our experiments in the Method section.In short, our active fluid is a flocking liquid assembled from colloidal Quincke rollers [3], which features bistable laminar flows in straight channels [13].Here, we confine the roller fluid in networks of channels.They form a honeycomb structure, which includes up to 3,200 trivalent nodes as illustrated in Fig. 1a.In their pioneering experiments Wioland et al. showed that active flows can strongly de- We can distinguish steady vortices in a finite fraction of the channels.(c) All the results reported in this article do not depend on the specific geometry of the nodes that connect the microfluidic channels.Top: node including a star-shape flow splitter.Bottom: Node with no splitter.In steady state only two of the three channels support laminar flows.The third (horizontal) channel hosts a steady vortex.(d) Distribution of the roller density in the channels.The density is homogeneous whatever the node geometry but features less fluctuations when including a star-shape splitter.(e) Distribution of the current Φe normalized by the most probable value Φ0 = 400 s −1 .We find a trimodal distribution for both node geometries.(f) The heatmap shows the local value of the current Φe supported by each edge of the channel network.The sign of Φe is arbitrarily defined using the bipartite nature of the honeycomb lattice.Φe is chosen to be positive when the flow goes from a given sublattice to the other.The solid line indicate the shale of the streamlines and the arrows the orientation of the flow.The data correspond to the same region of space showed in (b).(g) Classification of the seven possible vertices at each node.The color of the edges indicate the sign of the current.It is associated to a classical spin-1 variable Φe ∈ −1, 0, 1, (Φe = 0 means: no current on the edge e).We distinguish three different classes of vertices indexed by another spin variable σ ∈ −1, 0, 1 that distinguishes the handedness of the three-color patterns.σ is defined at the nodes of the network.σi = +1 when the set of colors at the node i is a positive permutation of {orange, blue, gray} indexed in the clockwise direction.σi = −1 for negative permutations, and σi = 0 when all three fluxes vanish, see also Methods.
pend on the specific geometry of their boundaries [23].
To establish the robustness of our findings, we therefore consider two different node geometries where the trivalent junctions either include a star-shape splitter or not, see Figs. 1c.We also systematically vary the aspect ratio = w/ of the roller steams in the channels, by varying their length from 120 µm to 200 µm, see Fig 1c .The channel width w is kept constant and equal to 200 µm, about a hundred times larger than the colloids' radius (2.4

µm).
We first observe a transient dynamics where the active fluid undergoes strong density and velocity fluctuations over system spanning scales, see Supplementary Video 1.It then reaches a steady state where the den- sity is homogeneous across the whole device, see Fig. 1a and 1d, and Supplementary Videos 2 and 3.By contrast, the mass-current distribution in the channels is heterogeneous, Fig. 1e.It is trimodal with two symmetric peaks associated to spontaneous flows at constant speed, and a higher peak revealing the presence of channels supporting no net current, see Fig. 1e and 1f.We find that the peaks of the distributions of the density and current (Φ e ) are narrower when using the splitter geometry shown in Fig. 1c.All the results reported below therefore correspond to this node geometry.
Unlike in disconnected straight channels, the active flows are never laminar throughout the whole sample: at least 40% of the channels host steady vortices and therefore support no net mass flux, see Fig. 1b, 1e and 1f.This series of observations stems from the geometrical frustration of active laminar flows.It can be understood as follows.As in most active fluids, the spontaneous flows of the Quincke rollers operate at constant speed.When the fluid density is homogeneous, the magnitude of the net fluxes |Φ e | measured along all the network edges (e) is peaked on a constant value Φ 0 (in our experiments Φ 0 = 400 s −1 ), Fig. 1e.The active flows are therefore frustrated at any node with an odd coordination.Flowing at constant speed and conserving mass are two incompatible constraints that cannot be simultaneously met in steady state.In the specific case of trivalent network, frustration is accommodated at each node by the emergence of vortices in either one, or three, edges where the flux vanishes (Φ e = 0), see Supplementary videos 4, 5 and 6.Irrespective of the channel aspect ratio, the situation where all three fluxes vanish is also possible but less likely to happen as seen in Supplementary Figure 1.
The geometry of the streamlines is therefore determined by the conflicting imperatives set by activity and mass conservation.The resulting local frustration defines a set of seven possible flow rules at the vertices which we classify in Fig. 1g.From a condensed matter perspective, they are akin to the spin-ice rules responsible for the ground state degeneracy of magnetic textures in frustrated magnets [8,9,20,21].More specifically, the six most probable vertices of Fig. 1g define a three-coloring model on the honeycomb lattice [25].This first analogy with spin-ice physics explains the vast degeneracy of the flow patterns found in our experiments.Active hydraulic flows are not deterministic.Repeating the same experiment in the same periodic geometry, we observe a plethora of disordered flow patterns.We illustrate them in Fig. 2 and show that they hardly feature any correlation.
The streamlines form nearly close-packed ensembles of self-avoiding loops.This geometry is a direct consequence of the spin-ice rules sketched in Fig 1g, where the seven vertices are the generators of the self-avoiding random walks on the honeycomb lattice.The closing of the streamlines into loops readily follows from mass conservation.We note that these patterns are consistent with the AFN model numerically discussed in [16,18] FIG. 3. Polymorphism of the streamlines.(a) Evolution of the streamline morphology with the channel aspect ratio.The color indicates the flow direction along the loops (Orange: clockwise, purple: counterclockwise).Top row: Experiments.The two set of streamline loops correspond to = 0.7 and = 1.1 respectively.Below = 0.8, the streamlines are crumpled and segregated.By contrast, above , they form nested and persistent loop patterns.The two heat maps are the topographic map of the same experiments.The local height field is defined as the sum of the winding numbers associated to all loops winding around a given point in space, see Methods.Below the segregation of the streamlines translate in a two-level topography.Conversely, above the map features high hills and deep valleys.Bottom row: Numerical simulations.Streamlines and topographic maps predicted by the low energy configurations of the Hamiltonian defined by Eq. ( 1) for JC = −2 and JC = 2.The qualitative agreement between the simulations and experiments confirm the relevance of our active-hydraulics laws.It also shows that JC is a good proxy for the channel aspect ratio.(b) Characterization of the streamline geometries: experiments, simulations and theory.All simulations are performed on a honeycomb graph having the same size and boundary conditions as in the experiments.First column: For all values of (experiments), and JC (simulations), the gyration radius Rg grows algebraically with the length of the streamline loops N : Rg ∼ N ν .The dashed lines are the two theoretical predictions detailed in SI based on loop O(n), and three-coloring models.Second column: Variations of ν with (experiments) and JC (simulations).Dashed-dotted line: exact prediction in the small JC and large JA limit (closely packed O(1) loop model).Dashed line: exact prediction for JC = 0 and JA 1 limit (closely packed O(2) loop model).Third column: Nesting level quantified by the maximal height difference on the topographic maps (average over ten realizations).The nesting level jumps at = (experiments), and JC = 0 (simulations) Fourth column: C(r) is the probability to find two nodes at a distance r within the same loop.When > the correlation C(r) decays much faster than when < .The trend is even clearer in the simulations, where the curves collapse on exponential and power law mastercurves.Dash-dotted line: Exact prediction in the limit JC 0, dashed line: exact prediction for JC = 0.
for degree-3-vertex graphs.However the loops' morphology strongly depends on the aspect ratio of the channels , Figs. 3a.Obviously, this polymorphism cannot be explained by the sole spin-ice rule of Fig. 1g which is agnostic to the channel lengths.To gain more insight on the origin of streamline patterns, we quantitatively characterize the loops' geometry.Whatever the channel aspect ratios, the loops are self-similar: their gyration radius R g grows algebraically with their length L, Fig. 3b.For small values, the loops are collapsed and segregated, by contrast, when is large, the streamlines are more persistent and form nested structures.We quantify these observations by plotting the exponent ν of the gyration radius and see that it undergoes a sharp increase when exceeds = 0.8, Fig. 3c.
The transition from segregated to nested loops occurs at the same value .To see this, we classically identify the orientated streamlines with the contour lines of the height fields h of rough landscapes [6].We show in Fig. 3a the corresponding topographic maps of h, see also Methods.We then quantify the level of nesting by the maximal height difference ∆h which increases sharply at , Fig. 3b.These observations prompt us to investigate deeper the orientational interactions between the streamlines.To do so, we measure the fraction of parallel configurations when two adjacent streamlines are separated by a channel supporting no net flux, Fig. 4b.To measure this quantity we note that it is given by 1 + σ 1 σ 2 /2 where the spin variables σ i measure the handedness of the vertices see Fig. 1g and Fig. 4a.
Fig. 4b shows that antiparallel contacts (σ 1 σ 2 = −1) prevail when < , whereas most of the contacts are parallel (σ 1 σ 2 = +1) when > .This central result indicates that active hydraulic flows are not only shaped by the spin-ice rules but also by short-range orientational interactions between adjacent streamlines.To elucidate their nature, we investigate the morphology of the vortical flows at the subchannel scale, Fig. 4a.We find a clear structural change at .Below , the channels with no net flux (Φ e = 0) host a single vortex.This vortex couples the adjacent channels as a gear would couple two circulators.The continuity of the flow field then favors antiparallel couplings between neighboring streamlines.This short range coupling is consistent with collapsed streamline loops including a number of hairpins, see Fig. 4e.
The situation where > is more subtle.We find that the zero-flux channels host two vortices, Fig. 4a.In most channels, they rotate in opposite directions, Fig. 4b.The continuity of the flow field therefore promotes parallel couplings, which explains the emergence of nested structure of the streamlines as clearly seen in Fig 4f .However in Fig. 4a we can also see configurations with two corotating vortices which promote antiparallel couplings.We now need to understand whether the dominance of counter-rotating vortices is specific to our experiments, or generic to polar active fluids.We address this ques-tion numerically.We model the active flows using Toner-Tu hydrodynamics of polar active matter [26,27], and solve these generic equations using a finite element solver, see Methods.We consider the simple geometry of an anisotropic cigar-shaped chamber with tangent boundary conditions, see Figs. 4c.A systematic investigation would go beyond the scope of this article, we therefore focus on a single set of material parameters consistent with earlier measurements in Quincke-roller fluids [28,29].As the channel aspect ratio increases, we observe a clear transition between one-vortex and two-vortex configurations.As a vortex hosts a +1 topological charge, when two vortices are present, an additional −1 topological defect must coexist with them to conserve the overall +1 topological charge [28].Remarkably, we find that the most likely configuration does not correspond to the coexistence of +1 and −1 charges, Fig 4d.By contrast, we observe a fractionalization of the −1 charge into two −1/2 singularities, Fig 4c .The fractional charges are stuck on the long channel walls and stabilize the counterrotation of the two +1 vortices.This numerical result is in excellent agreement with our experimental findings Fig 4a .We therefore conclude that the prevalence of parallel couplings between the macroscopic streamlines (Fig. 4b) originates from defect fractionalization in the flow field, and does not rely on the specifics of Quincke rollers.
We are now equipped to state the laws of active hydraulics and map them on a double spin model: (i) The currents Φ e on the edges e = {i, j} can only take three values: −1, 0, +1.It is a classical spin-1 variable.(ii) The stability of the laminar active flows in confined geometry penalizes the Φ e = 0 state.(iii) In steady state, mass conservation imposes the spin-ice constraint j Φ ij = 0 which frustrates (ii) at all nodes having an odd coordination number, Fig. 1g.(iv) Finally, the aspect ratio of the channels results in effective (anti)ferromagnetic interactions between adjacent streamlines.
To express these four elementary laws in more quantitative terms, we identify the streamlines with the equilibrium configurations of an effective Hamiltonian H that couples two spin variables: Φ e = ±1, 0 , and σ i = ±1, 0 which defines the handedness of the i th node (see Fig. 1g  and caption): The first term reflects the activity of the fluid.J A is a positive constant, which penalizes the states of vanishing flows (Law ii).The second term reflects the orientational interactions between adjacent streamlines (Law iv).J C is a proxy for the channel aspect ratio: J C < 0 corresponds to antiparallel couplings, i.e. to channels where < , conversely J C > 0 corresponds to parallel couplings, i.e. to channels where > , see Figs. 4a, 4e and  4f.We stress that J C only couples vertices connected by channels that support no net flux.We minimize H using a Monte Carlo worm algorithm [11] which enforces the spin-ice rule imposed by mass conservation (Law iii).We then compare the equilibrium states to our experimental measurements in Fig. 3b.The excellent agreement between the computed and measured streamline geometries confirms the predictive power of our active hydraulic laws.
We can now gain a deeper insight by exploiting quantitative analogies with a series of statistical field theories that were hitherto unrelated to any experimental system.Let us first discuss he limit of strong parallel interactions (J C 0), which is the most complex to analyze.This difficulty stems from the subexetensive degeneracy of the ground state.This observation was first made on the three coloring model introduced in Ref. [10] which maps on Eq. ( 1).The ground state configurations indeed correspond to fully nested patterns of streamlines, which we never find in our experiments and simulations.In fact, we probe metastable states for which no exact prediction is available.We however fully characterize them in SI, and show in Fig. 3b that our model correctly accounts for our experimental measurements in finite-size systems.
Mappings to field theories are more fruitful in the other two limits J C 0 and J C = 0.In these cases, our spin model allows exact predictions of the streamline geometry.In the limit of small aspect ratios, viz.J C 0, we can unambiguously define a net circulation around each face of the honeycomb lattice, see SI Fig. S9b.The handedness of the circulation around each face defines yet another classical spin variable Ω = ±1, which corresponds to 2|h| − 1, see Figs. 3a and 3b.Altogether, the constraints imposed by J A 0 and J C 0 translate into antiferromagnetic couplings between the adjacent Ω spins.As a result, in this regime, the streamlines are nothing else but the domain walls between regions of opposite magnetization in the ground states of a frustrated Ising antiferromagnet, see Fig. 3b and SI Fig. S9b.We use this analogy to provide exact expressions for the gyration radius, the nesting level and the probability C(r) to find two edges at a distance r in the same loop, Fig. 3b.The domain walls of the Ising antiferromagnet (AF) are obviously segregated and the nesting level is given by ∆h = 1 [10].The exact predictions of the gyration radius exponent and of the decay of C(r) are more complex to compute.They rely on the mapping of AF Ising models on the loop O(1) theory [17], which we recall in SI.Using this mapping, we find ν = 4/7 [21,24] and C(r) ∼ r −1/2 .These three exact predictions agree with our experimental findings, and numerical simulations, in the limit < (J C < 0), see Fig. 3b and SI.Finally, when J C = 0, the micro-scopic defect fractionalization does not affect the streamline orientations, and our model reduces to Baxter's three coloring model [25].It consists in coloring the edges of a honeycomb lattice as sketched in Fig. 1g.Baxter's model map on an exactly solvable statistical model of interacting loops known as the O(2) loop model, which we recall in SI as well.The loops of this model identify to the streamlines.Using this second powerful analogy we find that the gyration radius exponent should be ν 0 = 2/3 [21], and that C(r) ∼ r −1 .Remarkably these two exact predictions agree with the values measured in our experiments and simulations at the crossover between the parallel and antiparallel regimes Fig. 3c.Altogether the agreement between our experiments and theories establish the predictive power of our frustrated-spin model which accurately describes the random geometry of frustrated active-hydraulic flows.
As a final remark, we stress that the four local laws of active hydraulics apply broadly, beyond the specifics of periodic lattices and polar active matter.They should describe the flows of any form of active fluid animated by spontaneous laminar flows in complex networks, from cell tissues, to bacteria suspensions, to active gels and liquid crystals.We therefore expect our findings to provide a robust set of design rules for active microfluidic devices and offer new insights on the dynamics of groups of living cells and animals in heterogeneous environments and complex habitat [36][37][38].sodium salt).The microfluidic devices are made of two electrodes spaced by a 25-µm-thick double-sided tape.We pattern the bottom surface of the device with channel networks.They are made of a 2-µm-thick layer of insulating photoresist resin (Microposit S1818) patterned by means of standard UV lithography as explained in [1].In all our experiments the networks have the shape of a honeycomb structure, see Fig. S5a.The width of the channels is 200 µm and we vary their lengths from 120 µm to 220 µm.The overall size of the channel network is in the order of 1.5 × 1.5 cm 2 .
The electrodes are glass slides coated with ITO (indium tin oxide, Solems, ITOSOL30, thickness 80 nm).We start the experiments by filling the microfluidic chambers homogeneously with the colloidal solution, we then let the colloids sediment on the bottom electrode.The average packing fraction of the colloidal monolayer is approximately equal to 30% in all our experiments.We operate at this packing fraction to make sure that our active fluid remains in the polar-liquid state, far from the flocking-transition threshold (approximately equal to 0.1%) [3].We then apply a DC voltage of 110 V to trigger the Quincke instability and cause the colloids to roll at a constant speed v 0 ≈ 0.8 mm.s −1 on the bottom electrode.Each experiment is repeated from six to ten times.

B. Impact of the node geometry on the emergent flow patterns
We use two different designs for the nodes of the network.The first design is a simple junction between three straight channels, whereas the second design includes a flow splitter at each node, Fig. S5a.Our results do not qualitatively depend on the specific geometry of the nodes: the streamlines form self-avoiding loops, Fig. S5b, and the fraction of ferromagnetic couplings increases with the aspect ratio of the channels, Fig. S5c.However, including flow-splitters results in more homogenous fluid flows, Fig. S5d, and reduce the crossings between the streamlines, see Figs.S5b.In the main text, all the results correspond to the splitter geometry.

C. Measurement of the velocity fields
Once the active flow reaches a steady state, we image the whole network for 5 s with a Nikon AZ100 microscope.We record the videos with a Luxima LUX160 camera (Ximea) at a frame rate of 200 fps.To measure the velocity field v we use a conventional particle-imaging velocimetry (PIV).In practice, we use the PIVLAB MATLAB package [4].The PIV box size is 48 × 48 µm 2 , the PIV boxes overlap over a half of their size, Fig. S6a and S6b.Before constructing the streamlines, we average the velocity field over the width of each channel as shown in Fig. S6.We define the average velocity in each channel as a scalar quantity.We use the bipartite geometry of the honeycomb lattice to define its sign.Our sign convention is more easily explained using the sketch of Fig. S6c.Denoting by 'a' and 'b' the two sublattices, we choose to assign a + sign (blue color) to the velocity when the fluid flows from a 'a'-node to a 'b'-node and a − sign (orange color) when the fluid flows from a 'b'-node to a 'a'-node.

D. Density and current fields.
To measure the local current fields j(r) ≡ ρ(r)v(r), we measure the density and the velocity field in two different movies of the same experiment.For both measurements, we use a Hamamatsu ORCA-Quest qCMOS camera mounted on a Nikon AZ100 microscope with a 2× objective.As we cannot image the whole network at once, we perform multiple acquisitions and stitch our images, which is always possible in the steady state.To measure the density field, we perform epifluorescence imaging at 10 fps, and use the local fluorescence intensity as a proxy for the colloid density.Fig. S7a shows that, as expected, both quantities are proportional to one another.The calibration was performed on static images using a higher magnification and using a conventional particle detection algorithm (ImageJ).
We then use bright field imaging and record at a higher frame rate (200 fps) to perform our PIV analysis.The PIV box size is 24 × 24 µm 2 .We then average both density, ρ, and velocity fields v over time and we multiply them to reconstruct the local current field j see Figs S7b, S7c and S7d.Finally we measure the average current Φ ij along the edges of the network by averaging j over the channel joining the nodes i and j.

II. DATA ANALYSIS
In this section we detail how we analyse our data.We use the same methods and algorithms for our experimental and numerical data.

A. Construction of the streamlines
To construct the streamlines we first define a discrete current field Φ e on the edge e connecting the nodes i and j.Φ e can take three different values: ±1 in channels that support a net current and 0 in channel hosting vortices.In practice we construct the Φ e field as follows.We measure the average scalar product between the local velocity and the unit vector pointing in the direction of a channel.The average is performed over the channel area.When this quantity is larger than 0.5 we set Φ e = +1, when it is smaller than −0.5, we set Φ e = −1, and Φ e = 0 otherwise.In our simulations Φ = ±1, 0 by definition (see section III A).
Once the Φ e are defined we can unambiguously construct the oriented streamlines along the principal axis of the honeycomb lattice.We then use a depth-first search algorithm to detect and label each individual loop in the streamlines soup [5].Once the individual loops are identified, we can readily measure their gyration radius and the probability C(r) that two edges of the lattice separated by a distance r belong to the same loop (see section IV B).

B. Constructing a topographic map from oriented streamlines
We now explain how to quantify the nesting of the streamlines.To do so, we first convert the ensemble of oriented loops into a topographic map.We note F the ensemble of the hexagonal faces of the honeycomb network.Considering a configuration {Φ e } of the current field on the edges, we define a height field h f on the faces.
The loops are oriented, we can then define the winding number of each loop around a given face f , it is equal to ±1 when the face is lassoed by the loop and 0 otherwise.The height field of the topographic map is then defined at a face of the network as the sum of the winding numbers of all loops winding around this point.As a reference, the height is taken to be zero outside the network.
In practice we measure the height field h f recursively.The procedure is easily understood from the sketch of the algorithm shown in Fig. S8a.In short, starting from a face f associated with a height h we move to the neighbouring face f by crossing an edge e.If e hosts a spin pointing towards the right (resp.left) hand side, we assign the value h = h − 1 (resp.h = h + 1).The right and left directions are defined with respect to the vector connecting the centers of f and f .When the edge hosts no current (Φ e = 0), then h = h.The resulting height field does not depend on the way the network is explored.This procedure is very similar to the standard mapping of loop O(n) models to solid-on-solid models discussed e.g. in Ref. [6].
The reason why we introduce this mapping is twofold.Firstly, the height amplitude ∆h = max F h f − min F h f is one of our key observables.It quantifies the nesting of the streamlines and distinguish between the two phases observed both in the experiments and in the simulations, see see section IV B. Secondly, this mapping will also allow us to map the streamlines of the crumpled phase to the domain walls of the antiferromagnetic Ising model on the triangular lattice of faces F, see section V C.

III. THEORY: MODELING THE EMERGENT FLOW PATTERNS, LOOPS ON THE HONEYCOMB LATTICE A. Mapping active flows on a spin model
In this section, we introduce a spin model that effectively accounts for the geometry of the streamlines observed in our active hydraulic experiments.Constructing a predictive model from first principle interactions between the active colloids, or even from the Toner-Tu hydrodynamic theory would be a formidable challenge.Here, we instead build a phenomenological description from conservation laws and our robust experimental observations.

Honeycomb lattice geometry
In all that follows, we note E the set of edges of the honeycomb network, the channels, and V the set of vertices, see Fig. S9a.Each edge links two vertices, and each vertex in the bulk connects to three edges.In our numerical simulations, the networks are either be bounded, to make quantitative comparisons with our experiments (section IV), or periodic to investigate the asymptotic scalings of the streamline geometry (section V).It might be worth noting that the edges of the honeycomb lattice form a Kagome lattice.We call N the linear number of hexagonal faces (N = 3 in Fig. S9a), so that in a periodic system there are N 2 faces, 3N 2 edges and 2N 2 vertices.
Our model accounts for two constraints: (i) The active fluid current Φ e flowing along the edge e ∈ E can only take three values in each channel; and (ii) the mass is conserved at every vertex.

Active flows and Blume Capel spins
To model condition (i), we consider Blume-Capel spins Φ e ∈ {−1, 0, 1} [7] on the edges e = (i, j) ∈ E of the honeycomb lattice (which are the vertices of the dual Kagome structure).Φ e = 0 corresponds to the absence of net flows along the edge e: the channel hosts vortices.Φ e = ±1 corresponds to a finite current along e.The sign of Φ e indicates the direction of the flow as follows.The honeycomb lattice is bipartite: the vertices in V can be decomposed into two sublattices: 'V a ' and 'V b ', see Fig. S9a.The edges only link 'a'-vertices to 'b'-vertices.Given this observation, we assign Φ e = +1 to flows from an 'a'-vertex to a 'b'-vertex, and Φ e = −1 to flows from a 'b'-vertex to an 'a'-vertex.To distinguish between the three possible flow configurations, we associate a different color to each value of Φ e as illustrated in Fig. S9a: blue for Φ e = +1, orange for Φ e = −1, and gray Φ e = 0.

Mass conservation and spin-ice rule
In steady state, mass conservation imposes a hard constraint on the Blume-Capel spin configurations.The sum of the flows vanishes at each vertex, which translates into where ∂i is the subset of edges connected to the vertex i.In the bulk ∂i includes three spins, associated with seven configurations: {+1, −1, 0} = {blue, orange, gray} (and all possible permutations), or {0, 0, 0} = {gray, gray, gray} as illustrated in Fig. S9a.On the boundary of finite networks, there are two spins in ∂i, and they can only take the values {+1, −1}, {−1, +1} or {0, 0}.
We stress that the mass-conservation constraint of (S1) here corresponds to a spin-ice rule.As a matter of facts, in spin-ice models frustration arises from a similar constraint: the sum of the magnetic moments at each site of a crystalline spin lattice must vanish [8,9].We will show that the ice rule induces a strong degeneracy of the possible spin configurations in steady state.

Self-avoiding loops
As discussed in the main text, the spin-ice rule constrains the spin configurations, i.e. the streamlines of the active flow, to form oriented self-avoiding loops on E, Fig. S9a.This emergent geometry can be understood as follows.We start from an edge hosting a +1 spin (blue edge) pointing in the direction of the vertex i.The spin-ice rule imposes that this vertex is also connected to a single −1 spin (orange edge) pointing in the direction of a neighboring vertex.We can then iterate this reasoning to construct a walk on the honeycomb lattice.As each vertex can only connect to a +1 and a −1 edge the walks cannot intersect or cross themselves as these crossings would require vertices connected to three non-zero spins.Finally, the resulting walks necessarily form closed loops in virtue of mass conservation.
We note in passing that any ensemble of oriented self-avoiding loop on the honeycomb lattice can be tough as an ensemble of Blume-Capel spins on a Kagome lattice and obeying the spin ice rule.

Interactions between the loops
As discussed in the main text, the spin-ice rule is not sufficient to explain the the geometry of the active hydraulic flows observed in our experiments.To make progress, we need to account for the "parallel" and "antiparallel" flow couplings across non-flowing edges, Fig. S9b.We therefore add a minimal ingredient to our spin model.Fig. S9b would naturally suggest to include four-spin (anti)ferromagnetic interactions.Although possible this approach would lead to rather complex couplings.In order to simplify our model, we instead assign an additional spin variable σ i ∈ {+1, −1, 0} to each vertex i ∈ V visually defined in Fig. S9c.This spin variable describes the handedness of the Φ e spin configurations around a vertex.In short, we measure the values of the Φ e around the vertex i in the counter-clockwise direction.σ i = 1 when the the edge spins take the values (+1, −1, 0) (ie.blue, orange, gray) up to a circular permutation.σ i = −1 when the edge spins take the value (−1, +1, 0) (ie.orange, blue, gray).Finally σ i = 0 at vertices where the three edge spins are (0, 0, 0) or at the boundary of finite systems.Fig. S9c then readily tells us that (anti)parallel couplings between the streamlines translate into (anti)ferromagnetic interactions between the σ spins.We note that this construction is closely related to the color-dependent interactions of the 3-coloring model defined in Ref. [10].

Hamiltonian description of the stream lines
We can now define a phenomenological Hamiltonian model to account for the experimental flow patterns.We introduce the energy H of a spin configuration {Φ e }: where σ i and σ j are functions of {Φ e } and δ Φe,0 = 1 if Φ e = 0 and 0 otherwise.The two coupling constants J A and J C have a clear physical meaning.J A > 0 is a positive constant that promotes spontaneous flows, i.e that favors finite spin values Φ e = ±1.J C can be either positive, or negative, to account for the parallel, or antiparallel, couplings between streamlines separated by channels hosting vortices.

B. Numerical methods: Monte-Carlo worm algorithm
Our idea consists in describing the streamlines of the active hydraulic flows as the low energy states of H satisfying the spin-ice constraint (S1).To find them, we use a Monte-Carlo algorithm.We therefore define a statistical mechanics model given by the following Boltzmann weights at inverse temperature β, with the partition function OL(E) is the ensemble of oriented loops on the honeycomb lattice.Sampling the spin configurations on this ensemble guarantees that the spin-ice rule is satisfied.In our numerical simulations, we set J A = 1 without loss of generality and focus on the low temperature limit β 1.The main parameter of our model is then J C that regulates the coupling of the flows at the two ends of a non-flowing channel.To satisfy the spin-ice rule, we use a standard tool for lattice loop models, namely a worm Monte Carlo algorithm [11,12].In brief, the algorithm consists first in generating a worm, that is to say a loop of edges that can be flipped without violating the ice rule.The Monte Carlo move then corresponds to flipping the edge spins from ±1 to 0 or 0 to ±1 along this worm.In details, we use an adaptation of the short loop algorithm developed in Ref. [11] for the three-color model.We initialize the simulation with all edge spins equal to zero (Φ e = 0).This configuration trivially satisfies the ice rule.We then generate a new allowed configuration by constructing a worm.We explain the procedure to generate it in general, and not only for the initial configuration.Starting from a given spin configuration {Φ e } satisfying the ice rule, we choose an edge e at random, Fig. S10a.It is the first edge of the worm.If Φ e = ±1, we update it to Φ new e = 0, and if Φ e = 0 we update it at random to Φ new e = ±1, Fig. S10b.In both cases, this creates a pair of defects to the ice rule on the two ends of the edge e = (i, j): these defects are referred to as charges in spin ice models: Q i = e∈∂i Φ e = 0. To construct the second edge of the worm, the idea is then to propagate the charges on the lattice until they annihilate.To do so, we consider the site i and try to recover Q i = 0 by changing the spin of one of the two neighboring edges of i (different from e), Fig. S10c.When a single move is possible, we accept it and the second edge of the worm is defined.If there exists two possible moves, we choose one at random, which then defines the second edge of the worm.We iterate this procedure which results in the construction of a walk on the honeycomb lattice, Fig. S10d.When this walk visits a site k that has already been visited, we form a self-avoiding loop attached to a dangling end at k, Fig. S10e.We remove the dangling end by discarding all the spin flips used to construct the first steps of the random walk, Fig. S10e.This procedure is known as the short loop algorithm [11].
We can now evolve the Monte-Carlo dynamics by computing the energy from Eq. (S2).The spin-flips along the worm are then accepted according to the standard Metropolis rule [13] using the Boltzmann weight (S3) at inverse temperature β.This algorithm is ergodic (any configuration can be reached from the all-zero state and vice-versa) and satisfies detailed balance (Metropolis rule): it thus samples the ensemble of configurations compatible with the ice rule according to the Boltzmann probability distribution (S3).
We performed two types of simulations.(1) Small systems without periodic boundary conditions to compare our numerical results to our measurements.(2) Large systems with periodic boundary conditions in order to compute accurate scaling exponents.Simulations ( 1) are performed at system size N = 30, inverse temperature β = 2, flow constant J A = 1 and coupling constant J C = −2, −1, −0.5, 0, 0.5, 1, 2. For each set of parameters, we perform 50 independent simulations.We perform 10 6 worm iterations to thermalize the system, and then 10 6 additional iterations during which we compute the observables every 10 3 iterations.We compare our results with our experiments in Section IV.Simulations (2) are performed for system sizes ranging from N = 20 to N = 500, β = 1, J A = 1 and J C = −1, 0, 1.For each set of parameters, we perform 10 independent simulations, with 20000N 2 thermalization iterations and 20000N 2 more iterations during which the observables are computed every 200N 2 iterations.We report our results in Section V.

A. Qualitative features
We performed simulations of our spin model with 30 × 30 hexagons, which is close to the size of our microfluidic networks (see section III B).In the experiments, the control parameter is the aspect ratio of the channels.It controls the fraction of antiferromagnetic couplings.For < * = 0.8 the flow couplings are mostly antiparallel while for > * they are mostly parallel.The analog in the simulations is the coupling constant J C : J C < 0 favors antiparallel couplings, whereas J C > 0 promotes parallel couplings.
Fig. S11 confirms that the sign J C is a very good proxy for the role of the channel aspect ratio of the streamline geometry.We indeed find that when J c < 0 the spin loops are crumpled and collapsed, as the streamlines in the experiments where < .In addition, the equivalent topographic map typically includes only two to three values, again in excellent agreement with our experimental findings.In short channels, the streamline soup does not form a nested structure.When J c > 0, we find that the spin loops are more persistent and form nested structures, the height functions takes a range of values spanning a typical interval ∆h 8 for our system sizes.These geometrical feature again reproduce what we observe in our experiments where the streamlines are persistent and crumpled when > .In all cases our simulations also reflect the strong degeneration of the steady states observed in our active hydraulic networks.
This qualitative agreement begs us to quantitatively characterize the geometries of the two loops networks to validate the relevance of our analogy between frustrated Blume-Capel spins and active hydraulic flows

B. Quantitative comparison: geometry of the streamlines and spin loops
To be more quantitative we consider all the observables accessible to our experiments and proceed to a systematic comparison with our simulations, Fig. S12.
• Fraction of parallel contacts.An essential observable distinguishes the two active flow patterns: the fraction of parallel contacts across non-flowing edges f para = n para /(n para +n anti−para ) where n para is the number of edges hosting no net flow (Φ e = 0 in the simulations) and separating two streamlines flowing in the same direction (σ i σ j > 0 in the simulations), n anti−para counts the edges separating streamlines flowing in opposite directions.
In the experiments, f para varies from ∼ 0.1 for < * to ∼ 0.6 for < * , see Fig. S12e.i.Anti-parallel contacts dominate in the phase where the loops are crumpled while they are subdominant in the phase where the loops are more persistent and nested.This distinction is even more marked in the simulations: f para is very close to 0 for J c < 0 while it reaches 0.9 when J c > 0 (Fig. S12e.ii).
• Gyration radius.To distinguish between the crumpled loops of the first regime and the more swollen loops of the second regime, we compute the gyration radius R g of each loop, where the sum runs over the L edges e, at positions X e that make the loop.This is a typical observable for polymers [14] or spin-ice models [12].We consistently find that the loops have a self-similar geometry, the average gyration radius scales as a power-law of the loop size, In Fig. S12a.i and a.ii, we plot the gyration radius of the loops.Our experiments and simulations are consistent.Loops in the parallel regime have a gyration radius that grows faster than those in the anti-parallel regime.The curves collapse on two master curves in the two distinct regimes.We plot the variations of the exponent ν in Figs.S12d.i and d.ii as a function of and J C respectively.The transition is sharp both in the simulations and the experiments.The anti-parallel regime corresponds to an exponent ν that is consistent with 4/7 0.571.This exponent correspond to the gyration radius of the domain walls in the anti-ferromagnetic Ising model, see subsection V B. In the parallel regime we find a higher exponent, around 0.75 − 0.8 in the experiments and 0.82 in the simulations.We will argue below that there exist no indication for universality in this regime.
• Loop correlation.A typical observable that distinguishes loop models [12,15] is the probability C(r) that two edges separated by a distance r belong to the same loop.C(r) decays as a power-law in the critical phases of loop models.We plot this quantity in Fig. S12b.i and b.ii.We see that both in the experiments and the simulations, C(r) decays faster in the anti-parallel phase than in the parallel phase.The system size does not allow us to conclude on a power-law decay and on its exponent.However, the data in the anti-parallel regime are consistent with C(r) ∼ r −1/2 predicted for the anti-ferromagnetic Ising model (see section V B).Moreover, the data in the parallel regime may be consistent with C(r) ∼ r −1 (see section V D).
• Average height amplitude.Finally, we quantify the nesting level of the loops and plot the height amplitude ∆h of the patterns in Figs.S12f.i and f.ii.We see that ∆h undergoes a sharp jump at the point where the fraction of parallel contact jumps.This observation holds both in our experiments and simulations, as anticipated from the patterns plotted in Fig. S11.
This ensemble of results confirm that our minimal spin model capture all the salient features of the active hydraulic flows measured in our experiments.The size currently accessible to our experiments does not allow to predict whether the transition between the crumpled and the nested streamline regimes is a genuine dynamical phase transition or a smooth crossover.Answering this question goes beyond the scope of this article.However, we can take advantage of our numerical simulations to predict the asymptotic scaling laws that characterize the self-similar shape of the interacting self-avoiding loops.

V. ASYMPTOTIC LIMITS OF THE SPIN MODEL AND STREAMLINE GEOMETRY
In the previous section, we introduced an equilibrium spin model (Eqs.(S2)-(S4)) to account for the flow patterns observed in the experiments.We found a good agreement between experiments and simulations in both regimes of antiparallel and parallel flow couplings, whcih respectively correspond to small and large channel aspect ratio.To gain a more quantitative insight into the morphology of the streamlines, we investigate numerically the two asymptotic limits: J C → −∞ and J C → +∞ (in the so called fully-packed limit where J A → ∞).We use theoretical arguments and extensive numerical simulations to predict the scaling laws of the gyration radius ν and of the decay of the correlation C(r).Throughout this section, all simulations are performed with periodic boundary conditions.This technical section is organized as follows: Section V A is a reminder of the loop O(n) models and of their mapping on spin models.Section V B focuses on the antiparallel regime (J C → −∞).We show that this regime maps on the antiferromagnetic Ising model on a triangular lattice, or equivalently the loop O(1) model.This mapping allows us to predict the gyration radius exponent (ν = 7/4) and the algebraic decays C(r) ∼ r −1/2 .In Section V C, we investigate the case where there is no orientational coupling between the streamlines (J C = 0).While this situation cannot be realized in our experiments, it provides an insight into what happens in the crossover between the two asymptotic regimes (J C 0 and J C 0).When J C = 0 our model maps on a three-color model, or equivalently to a loop O(2) model.In this regime we predict a gyration radius exponent ν = 2/3 and a correlation decaying as C(r) ∼ r −1 consitent with our experimental findings at intermediate aspect ratio.Finally, in Section V D, we discuss the parallel regime (J C 0).This regime is not as simple as the two others as the ground state of our model features a sub-extensive degeneracy.The simulations reach metastable states which we characterize numerically.

A. Mapping spin statistics on loop O(n) models
The spin configurations constrained by the spin-ice rule correspond to oriented loops on the honeycomb lattice.Loop models have a long history in equilibrium statistical mechanics.A classic class of loop models is the loop O(n) models [15,16].We briefly review some key results as this class of model offer powerful tools to explain the asymptotic limits of our spin model, and the corresponding active-flow geometry.
Definition of the loop O(n) model.Let us consider the set of edges E of a finite lattice, in our case the honeycomb lattice.We define L(E), the set of all loop configurations on E. A loop configuration ω is a subset of E for which all the vertices of the underlying graph have an even connectivity: in our case each vertex is connected to either 0 or 2 edges.We denote #edges(ω) the number of edges in a loop configuration ω, and #loops(ω) the number of loops.By definition, the probability of ω in the loop O(n) model is given by with the partition function where x is the fugacity of the number of edges and n is the fugacity of the number of loops.The limit x → ∞ corresponds to the fully-packed regime: the number of edges in the loop configuration is maximized.
Three important cases: n = 0, n = 1 and n = 2.It is known that the loop O(n) model [15] maps on canonical statistical mechanics models.The limit n → 0 corresponds to the celebrated De Gennes mapping on the statistics of self-avoiding random walks.The case n = 1 maps on the Ising model on the triangular lattice of faces F, with x = e −2β where β is the inverse temperature of the Ising model.We will detail further and use this equivalence in Section V C. Finally, n = 2 at x = ∞ is equivalent to a three-coloring model which we discuss in Section V B to investigate the J C = 0 phase.

Relation between spin and loop O(n) models.
We recall that for integer n the loop O(n) model is closely related to the spin O(n) model [15,17].Considering n-dimensional spins s i = (s 1 i , . . ., s n i ), with k s k i = n, on the vertices i ∈ V of the graph, the partition function of the spin O(n) model reduces to Writing the weight function as w(s i • s j ) = 1 + xs i • s j , one can show that Z spin (x, n) = Z loop (x, n) [15,17].
Criticality and Schramm-Loewner evolution.One important question for the loop O(n) model, for a given loop fugacity x, is whether the loops have a characteristic size or not.In other words, does the probability of finding a loop of length L, P loop (L), decay exponentially or algebraically?This question was answered by Nienhuis in [6].He showed that for n ∈ [0, 2] there exists a critical point at x c (n) = 2 + (2 − n) 1/2 −1/2 : when x < x c P loop (L) decays exponentially (sub-critical phase), whereas for x > x c and x = x c P loop (L) decays as a power-law (critical phase).When x > x c , the loops are non only scale-invariant but conformaly invariant.It was then conjectured in [18] than in this critical phase, the scaling of the loop size is given a Schramm-Loewner evolution SLE(κ) [19] with parameter κ ∈ [4,8] given by n = −2 cos(4π/κ). (S10) In other words κ = 8 for n = 0, κ = 6 for n = 1, and κ = 4 for n = 2.The technical aspects of these results go beyond the scope of our article.We simply recall that a Schramm-Loewner curve of parameter κ [19] corresponds to a stochastic evolution in the half-plane having a fractal dimension d f and that features conformal invariance, In particular d f = 7/4 for κ = 6 (n = 1); d f = 3/2 for κ = 4 (n = 2).The correspondence between loop O(n) models and Schramm-Loewner evolutions (or more precisely conformal loop ensembles) is the topic of numerous works in the mathematical litterature but rigorous results remain scarce [15,16].
In the following two sections, we show that the phases J C < 0 and J C = 0 of our model correspond to fully packed loop O(1) and O(2) models.

B. Antiparallel flow couplings: antiferromagnetic Ising model on the triangular lattice
In this subsection, we investigate our model (Eqs.(S2)-(S4)) in the limit J C 0, i.e in the limit where all adjacent streamlines flow along opposite directions.This limit is directly relevant for experiments where the channel aspect ratio is small.

Mapping of our model on the antiferromagnetic Ising model in the limit JC → −∞
We focus on the ground states, β → ∞, of our spin model when J C < 0, with J A = 1.We show that they map on the ground states of the antiferromagnetic Ising model on the triangular lattice.
We first note that the constraint of maximal flow βJ A → ∞ implies that all vertices have their three edges associated with spins {+1, −1, 0}, there are no {0, 0, 0} vertices.This is readily seen by decomposing the sum where ∂i is the set of the three edges adjacent to i.This sum is maximized if all the internal sums are equal to 2, hence all vertices must have a {+1, −1, 0} structure.
Second, the limit βJ C → −∞ implies σ i σ j = −1 whenever a non-flowing edge (Φ e = 0) connects the sites i and j.Adjacent stream lines flow in opposite directions.Constructing the topographic map of subsection II B, we then readily find that the height field can only take the values h = 0 and h = +1, or h = 0 and h = −1, see Fig. S13.We can make yet again another spin analogy.For each face f ∈ F, we associate the two possible height values to Ising spins S f = +1 when h = 0 and S f = −1 when h = +1.
Finally we can show that the limit βJ C → −∞ implies antiferromagnetic couplings between the S f spins.This is done by noticing that if f and f are the two faces in contact with an edge e, then we can relate the edge-spin value to the face-spin values as Φ2 e = (1 − S f S f )/2.This equality follows from the rule at the origin of the construction of the topographic map.Using this simple relation we can recast the term that promote active flows on each edge as: which is nothing else but the energy of an antiferromagnetic Ising model where J Ising = J A /2 > 0. In the limit βJ A → +∞ the geometry of the streamlines therefore correspond to geometry of the domain walls in the ground state of an antiferromagnetic Ising model on the triangular lattice.This mapping readily explains the extensive degeneracy of the ground states which is a well established result for antiferromagnets on the triangular lattice, see e.g.[20].
In the next section we take advantage on this last frustrated-magnetism analogy to make exact predictions about the streamline geometries of active flows in trivalent hydraulic networks.for system sizes N = 20 to 1000, excluding loops which wind around the periodic boundaries.The exponents of the power laws found in the two sets of simulations are respectively 0.572 and 0.576 in agreement with the exact prediction ν = 4/7 ≈ 0.573.In the inset, the average lengths of the loops that wind around the system is plotted as a function of the system size.The growth is in agreement with N 7/4 (black line) which is the fractal dimension d f = 1/ν.(b) and (e) show the decay of the probability C(r) that two edges at a distance r belong to the same loop.The two sets of simulations are consistent with an algebraic decay as r −1/2 .(c) and (f) Distribution of loop sizes L (excluding winding loops) for the two types of simulations.The prediction for the exponent is τ = 2 + 1/7 ≈ 2.143.For all three observables, the agreement with theoretical predictions is more accurate for the Ising simulations.The statistics and the equilibration of the system are both better in this simulation set.
We identify the domain walls as the edges between two sites with opposite spins.The loops are found using a standard depth-first-search algorithm [5].The statistics of the loop lengths and the correlations between edges in the same loop are easily obtained.The gyration radius of a loop in periodic space is computed as follow.We consider a loop of L edges centered at positions (x 1 , y 1 ), . . ., (x L , y L ).We define x 1 = 0 and x i = x i−1 + per(x i − x i−1 ) for i = 2, . . ., L where 'per' is the shortest distance in periodic space.In other words, we redefine the coordinates iteratively starting from the first site of the loop.We proceed similarly for the y coordinates.The gyration radius R g of the loop is then defined as

C. Geometry of streamlines with no orientational couplings: 3-coloring problem and fully-packed loop O(2) model
We now turn to the limit J C = 0 of our active-hydraulics model.When J C = 0 the streamline loops are selfavoiding but their orientations are uncoupled.This situation is not directly observed in our experiments where the flow couplings are either parallel or antiparallel.However, our finite-size simulations indicate that in the crossover between the two regimes, when the fractions of parallel and antiparallel couplings are equal, our Φ e spin correctly accounts for the stream line geometry when taking J C = 0.
In order to gain a deeper insight into this limit, in the limit of asymptotically large system sizes, we proceed to another analogy with frustrated statistical mechanics.We map the statistical mechanics of our model defined by Eq. (S2) on the 3-coloring problem, and the fully-packed O(2) model.We then take advantage of exact predictions from conformal field theory and compare them with our numerical simulations and experimental measurements.

Active hydraulics, 3-coloring model and loop O(2) model: exact results
We proceed as in the previous section, we focus on the maximally flowing regime: βJ A → ∞ in Eq. (S2).In this limit the vertices of our Φ e -spin model are coupled to edges with spin configurations of the type {+1, −1, 0} = {blue, orange, gray}.Since there is no additional constraint constraint (J C = 0), our model then exactly corresponds to a 3-coloring model defined on the edges of the hexagonal lattice the edges, see Fig. S15a.This classical statistical mechanics model was first introduced by Baxter in Ref. [25].It was later shown this 3-coloring model is equivalent to the fully packed loop O(2) model [21].We recall that n can be interpreted as the fugacity of the loops in a loop model.The value n = 2 can then be intuitively understood as there are two ways to orient each loop whose conformation are unaffected by their orientation when J C = 0.This situation contrasts with the limit J C 0 where the orientation of a given loop fully determined the orientation of all the others, thereby leading to a fugacity n = 1.
The 3-coloring on the honeycomb lattice, and the fully-packed loop O(2) model have been extensively studied in the mathematical physics community, and a number of exact results are available.For instance, for the 3-coloring model, the degeneracy of the number of configurations Ω is known to scale extensively with the number N v of vertices: ln Ω = N v ln W with W ≈ 1.20872 [25].We also recall that the fully-packed loop O(2) model is commonly believed to follow the same statistics as the SLE(κ = 4) process [18].We build on this knowledge to make quantitative predictions on the stream line geometry and test these exact predictions against numerical simulations of our Φ e -spin model and experiments.
• Gyration radius.The fractal dimension of the loops in the fully-packed loop O(2) model is d f = 3/2 [21]: the gyration radius scales as R g ∼ L ν with ν = 1/d f = 2/3.This prediction in very good agreement with the numerical simulations of our spin model at J C = 0, see Fig. S16a, and is consistent with the value of ν measured in the crossover regime in our experiments.
• Loop correlation.The fully packed loop O(2) model, which has central charge c = 2 features algebraic correlations.The probability C(r) that two edges at distance r from one another belong to the same loop scales as C(r) ∼ r −1 [24].This result agrees with our numerical simulations, see Fig. S16b.
• Distribution of loop sizes.The probability P site (L) that the loop passing though a given vertex has size L scales as P site (L) ∼ L −(τ −1) with τ = 7/3 [21,24].The probability P loop (L) that a loop has size L scales as P loop (L) ∼ L −τ .These two predictions again agree with our numerical simulations as seen in Fig. S16c.theoretic terminology, we instead perform numerical simulations and address the similarities and differences with the other two phases (J C < 0 and J C = 0).We investigate the fully-packed loop configurations (2/3 of the edges are flowing) in the limit J C 0 where all flow couplings are parallel (σ i σ j = +1 accross edges hosting a zero spin Φ e ).The only solutions correspond to lines winding around the periodic box, and parallel to one another as illustrated in Fig. S17a.In a lattice of N × N hexagons, the number Ω of such configurations is of the order of 2 N [10].This can be seen by starting with a set of streamlines, or Φ e spins pointing along the same direction, say from left to right.To deform the loops, we can start from the leftmost edge, see Fig. S17 and flip the direction of the streamline, the second edge must host a current flowing to the right hand side, the streamline can then bifurcate upwards or downwards and so forth until the streamline winds around the box.This procedure implies that the number of configurations scales as the number of directed random walks which yields a 2 N scaling.This reasoning shows that the ground state degeneracy (more accurately the ground state complexity) is sub-extensive in the number of vertices N v = 3N 2 (or number of edges) since ln Ω ∼ √ N v ln 2. This behavior contrasts with the ln Ω ∼ N v ln W scaling found for the 3-coloring model when J C = 0, and with the algebraic scaling of the frustrated antiferromagnetic Ising model.
We stress that in a non-periodic system, the lowest energy configurations correspond to fully nested loops, whose degeneracy is also sub-extensive, following the same reasoning.However, both in periodic and finite systems these configurations are very hard to reach starting from random initial conditions.At low temperature, the worm Monte Carlo algorithm reaches highly stable local energy minima, see Fig. S17b.The same phenomenon happens in the experiments in which we never observed a set of perfectly nested streamlines configurations, see Fig. S12e.iwhere f para ≈ 0.6 is still much lower than 1.A full study of this potentially glassy phase is beyond the objectives of this article, but we provide numerical observations which shed some light on the geometry of the streamlines, i.e of the ground state of the Φ e -spin model.

Streamline geometry and non-critical phase for parallel couplings
We investigate again the same geometrical quantities numerically.
• Gyration radius.If we were probing ground state configurations, the gyration radius exponent would be trivially given by: ν = 1.In fact, the geometry which we observe is more complex.In the metastable states that we probe, Fig. S18a, we can distinguish two asymptotic behaviors: for short loops (L < 100), the exponent is ν 1 0.85 while for long loops (L > 100) the gyration radius has a scaling ν 2 0.6.The impossibility to describe the loops with a single exponent reflect that we are not dealing with a realization of a critical loop model.In the experiments, we probe loops up to size 300.This explains why we find an effective exponent ν exp 0.75 − 0.8, see Fig. S12d.i.
• Loop correlation.The probability C(r) that two sites as a distance r belong to the same loop is consistent with a r −1 decay (Fig. S18b).This scaling is identical to the simulations where J C = 0, see Section V C. But this scaling is a poor marker of different universality classes.The same scaling would indeed be observed in the ground state where all the streamlines would be parallel.Indeed, considering a given edge in a ground state configuration, Fig. S18a, there are 2 edges at a distance r that are in the same loop, out of a total number of edges that scales as r ( perimeter of a circle of radius r scales).In the ground states we therefore also find C(r) ∼ r −1 .Our numerical simulations tell us that the same scaling holds for the metastable states that we probe.
In conclusion, the metastable states of the parallel regime (J C > 0) are characterized by a gyration radius exponent ν 1 0.85 for short loops and a loop correlation C(r) that seems to decay as r −1 .Both observations are consistent with the experiments (Fig. S12 c.i and d.i).

VI. FINITE ELEMENT SIMULATIONS
Our goal is to determine whether the defect fractionalization seen in the flow field at the subchannel scale is specific to colloidal-roller fluids or generic to polar active matter.To answer this question we resort to numerical resolutions of Toner-Tu Hydrodynamics.Toner-Tu equations, Eqs.(S15) and (S16) are the equivalent of the Navier-Stokes for polar active fluids, flocks of self-propelled bodies [26,27].In their simplest form they reduce to [28]: where ρ(r, t) and v(r, t) are the local density and velocity fields of the polar active liquid.We choose the hydrodynamic constants to match the experimental values estimated in Refs [2,28,29]: λ = 0.7, σ = 5 mm 2 s −2 , D ≈ 10 −3 mm 2 s −1 , α 0 = 10 2 s −1 , ρ c = 3 × 10 −3 , β = 10 mm −2 s (we opt for a definition where ρ is nondimensional as in [28]).The physical meaning of each term has already been thoroughly discussed in the litterature, see e.g.[27].To solve this set of partial differential equations, we use the Finite Element Method described in Ref. [28].In practice, we use the open source software package FENICS [30].
To investigate the structure of the vortical flows in channels of different anisotropy, we solve the equations in cigarshaped boxes.The shape consists of a rectangular body with two semicircular segments attached to opposite sides.The rectangular body has a width of W and a length of L. The two semicircles have a diameter equal to W .We use the Babuska penalty method, (see [31]), to set tangent boundary conditions for the velocity field: v • n ≈ 0 at each boundary, where n is the normal vector to the mesh boundary.There is no boundary condition for the density field ρ, but we checked that mass is globally conserved.In all simulations, the density is initialized as a constant ρ 0 = 0.3, and the velocity vectors are initialized with a random orientation and a norm close to α/β.The time step between two time increments δt is chosen to be small compared to the typical relaxation timescale of the fast-speed mode τ = α(ρ 0 ) ≈ 10 −1 s.In practice, we take δt ≈ 10 −3 s.The computational mesh consists of ∼ 2000 triangular cells.In our simulations L varies from 0 to 200 µm while W is kept constant and equal to 200 µm.We interpolate the density and velocity fields using second-order polynomials on Lagrange finite-element cells (see [30]).

FIG. 1 .
FIG. 1. Active hydraulics in trivalent networks (a) Macroscopic picture of a trivalent microfluidic network.It encloses a two-dimensional active fluid made of flocking Quincke rollers.(b) Close-up view.We can distinguish steady vortices in a finite fraction of the channels.(c) All the results reported in this article do not depend on the specific geometry of the nodes that connect the microfluidic channels.Top: node including a star-shape flow splitter.Bottom: Node with no splitter.In steady state only two of the three channels support laminar flows.The third (horizontal) channel hosts a steady vortex.(d) Distribution of the roller density in the channels.The density is homogeneous whatever the node geometry but features less fluctuations when including a star-shape splitter.(e) Distribution of the current Φe normalized by the most probable value Φ0 = 400 s −1 .We find a trimodal distribution for both node geometries.(f) The heatmap shows the local value of the current Φe supported by each edge of the channel network.The sign of Φe is arbitrarily defined using the bipartite nature of the honeycomb lattice.Φe is chosen to be positive when the flow goes from a given sublattice to the other.The solid line indicate the shale of the streamlines and the arrows the orientation of the flow.The data correspond to the same region of space showed in (b).(g) Classification of the seven possible vertices at each node.The color of the edges indicate the sign of the current.It is associated to a classical spin-1 variable Φe ∈ −1, 0, 1, (Φe = 0 means: no current on the edge e).We distinguish three different classes of vertices indexed by another spin variable σ ∈ −1, 0, 1 that distinguishes the handedness of the three-color patterns.σ is defined at the nodes of the network.σi = +1 when the set of colors at the node i is a positive permutation of {orange, blue, gray} indexed in the clockwise direction.σi = −1 for negative permutations, and σi = 0 when all three fluxes vanish, see also Methods.

FIG. 2 .= 1 −
FIG. 2.Degeneracy of the streamline patterns.Two realizations of the same experiment ( = 0.7) lead to two different streamline patterns with markedly different geometries.See Fig.S7in SI for more realizations.In all realizations the streamlines form a soup of self-avoiding loops colored according to their length (the darker the longer).We stress that the overlap distribution computed over ten independent experiments indicates very little correlations from one realization to another.The overlap function O αβ quantifies the difference between the flow patterns in the two realizations α and β of the same experiment: O αβ = 1 − |Φα(e) − Φ β (e)| e, where Φ(e) = ±1, 0 is the sign of the flux in the edge e of the honeycomb graph, see Methods.The typical value of the overlap is 0.2, which is close to the value expected for uncorrelated random current fields.In the limit case where the Φα(e) are uncorrelated, we can perform a crude approximation.Ignoring spatial correlations, we define the probability 1 − 2p that Φ(e) = 0 and assume that Φ(e) = ±1 have the same probability p.The average overlap between two random configurations would then be O αβ e = 1 − 4p(1 − p).Estimating p from the fraction of edges supporting no net current, we find p ∼ 0.25 − 0.3, which yields O αβ e = 0.25 − 0.16.

FIG. 4 .
FIG. 4. Orientational interactions between adjacent streamlines and topological defect fractionalization.(a)Closeup on the flow field in two adjacent streamlines separated by a channel that supports no net flux.Below the central channel hosts a single vortex.Above the central channel hosts two vortices.They can either rotate in the same or opposite directions.When they rotate in the same direction, the flow field includes a −1 topological charge separating the two +1 charges at the center of the vortices.Conversely, when the two vortices rotate in opposite directions the −1 charge is fractionalized into two −1/2 defects bound to the channel walls.(b) Fraction of parallel contacts across each channel supporting no net flow plotted against the aspect ratio . is defined in Fig.1for the experiments, and is given by the ratio between the length and the width of the anisotropic simulation boxes.The fraction of parallel contacts is equal to1  2 (1 + σ1σ2), where σ1σ2 = 1 when the two vertical flows are along the same direction and −1 otherwise.The flow couplings are mostly antiparallel below and mostly parallel above .(c) Numerical resolution of the Toner-Tu equations in cigar-shaped boxes, see SI for details about the numerical method and material parameters.We denote w the width of the box, the length of the straight line joining the half-circles, the aspect ratio of the boxes is defined by = /w.In agreement with what we observe in b, increasing the aspect ratio of the anisotropic box result in a transition form a flow field with a single +1 topological charge to flow fields including two vortices, which can be either separated by a single +1 defect or two fractionalized −1/2 defects bound to the walls to conserve the overall topological charge.(d) Repeating the same simulation with random initial conditions we find that the steady states where the rightmost (v2) and leftmost (v1) flows are parallel prevail in boxes with large aspect ratios.In other words the transition from antiparallel to parallel contacts in our experiments are not specific to colloidal rollers but generic to flocking matter (i.e Toner-Tu fluids).(e) Close-up on the experimental streamlines measured in a channel where < .Close-up on the experimental streamlines measured in a channel where > .Orange : streamline loops flowing in the clockwise direction, purple : streamline loops flowing in the counterclockwise direction.The green arrow indicates the local orientation of the flows.The streamlines are crumpled and segregated.The antiparallel couplings between adjacent streamlines promote this geometry.(f) Same as in (e) for > .The streamlines are persistent and nested.The parallel couplings between adjacent streamlines promote this geometry.
FIG. S5.(a) Zoom on the two junctions geometries.Top: simple trivalent junction.Bottom: trivalent junction with a splitter.(b) Top: streamlines in a networks with simple trivalent nodes.= 0.7.Bottom: streamlines in a networks with splitters.= 0.7.Streamlines form self avoiding loops, but splitters forbid crossings between streamlines.(c) Fraction of parallel couplings for simple nodes (top) and for nodes with splitters (bottom).(d) Distribution of the average density in channels.Top : simple nodes, = 0.84.Bottom: nodes with splitters, = 0.84.(e) Distribution of the average current in channels.Top : simple nodes, = 0.84.Bottom: nodes with splitters, = 0.84.
FIG. S6.(a) Velocity field.The averaging windows are sketched in red.(b) Sketch of the averaged velocity in each averaging window.(c) Scalar value of the velocity.It is +1 (resp.−1) when the fluid flows on average from a 'a'-node to a 'b'-node (resp.from a 'b'-node node to a 'a'-node).
FIG. S7.(a) Fluorescence intensity with respect to the packing fraction ρ.(b) Colormap of the packing fraction field.(c) Velocity field.(d) Current field.

FIG. S8 .
FIG. S8.Mapping on a topographic map.(a) Configuration of our model (see Fig. S9).(b) Rules for the topographic map.If there is no arrow, the two faces have the same height.If there is an arrow pointing to the right, when moving from f to f , then h f = h f − 1.If there is an arrow to the left, then h f = h f + 1. (c) Topographic map corresponding to (a).In this configuration, ∆h = max h − min h = 2.
FIG. S10.Worm algorithm.(a) We select an edge in the initial configuration.(b) Flipping the spin (here 0 → −1) leads to the creation of a pair of defects (white circles).(c) We move one of the defects by flipping another edge spin.(d-e)We iterate the procedure until we come back to a site we already visited.(f) The beginning of the path is discarded.We use the Metropolis rule on the move that goes from configuration (a) to configuration (e).

FIG. S11 .
FIG. S11.Degeneracy of the flow patterns.Counter-clockwise loops are shown in black, and clockwise loops in green.The colors of the faces correspond to their heights (see subsection II B and Fig. S9c).Top: antiparallel regime with experimental configurations at = 0.7 and numerical configurations at JC = −2.Middle: numerical configurations with no flow coupling (JC = 0).Bottom: parallel regime with experiments at = 1.3 and numerics at JC = 2.

FIG. S12 .
FIG. S12.Quantitative comparison between numerical and experimental results.(a.i) and (a.ii)Average gyration radius of a loop as a function of its length for experimental aspect ratio ∈ [0.55, 1.44] (in (a.i)) and numerical coupling constant Jc ∈ [−2, 2] (in (a.ii)).(b.i) and (b.ii) Experimental and numerical probability C(r) that two sites at a distance r belong to the same loop.(c.i) and (c.ii) Same curves in log-log scale.The dashed and dash-dotted lines respectively correspond to exponents −1 and −1/2.(d.i) and (d.ii) Exponent ν of the gyration radius, fitted from (a.i) and (a.ii).The error bars are obtained from the least-square fit algorithm with input uncertaintly given by the spread in each bin of (a).The horizontal dash-dotted line corresponds to the prediction of the antiferromagnetic Ising model (4/7, see subsection V B). (e.i) and (e.ii) Fraction fpara or parallel streamlines across non-flowing edges.(f.i) and (f.ii) Height amplitude ∆h (see Fig.S9c).

) ν = 1 / 2
corresponds to Brownian motion, ν = 1 corresponds to straight lines or circles and ν = 3/4 corresponds to self-avoiding walks in 2D.Note that ν = 1/d f where d f is the fractal dimension of the model.
FIG. S13.Antiparallel regime (JC < 0): Mapping our active flow model on the antiferromagnetic Ising model.(a) Example of fully-packed configuration of our spin model (Eq.(S2)) with antiparallel couplings only (σiσj = −1 on all zero edges).The spin arrows reflect the local direction of the active flow in the fluidic netowrk (b) Mapping of the spin configuration of (a) on a topographic map (see Section II B).We note that the height field can only take two different values, which we associate to two Ising spins + and − to them.In this limit, the streamlines of the streamlines of the active hydraulic flows map on the domain walls of an antiferromagnetic Ising model.(c) Another perspective on our statistical mechanics problem consists in looking at the shape of the closed streamlines irrespective of their orientation (dark lines).They define a realization of a loop O(1) model.
FIG. S14.Antiparallel regime (JC < 0): numerical results.(a-c) Monte-Carlo simulations of the Ising model.(d-f) Monte-Carlo simulations of our spin model with JC = −2.(a) and (d) show the gyration radius Rg as a function of loop size Lfor system sizes N = 20 to 1000, excluding loops which wind around the periodic boundaries.The exponents of the power laws found in the two sets of simulations are respectively 0.572 and 0.576 in agreement with the exact prediction ν = 4/7 ≈ 0.573.In the inset, the average lengths of the loops that wind around the system is plotted as a function of the system size.The growth is in agreement with N 7/4 (black line) which is the fractal dimension d f = 1/ν.(b) and (e) show the decay of the probability C(r) that two edges at a distance r belong to the same loop.The two sets of simulations are consistent with an algebraic decay as r −1/2 .(c) and (f) Distribution of loop sizes L (excluding winding loops) for the two types of simulations.The prediction for the exponent is τ = 2 + 1/7 ≈ 2.143.For all three observables, the agreement with theoretical predictions is more accurate for the Ising simulations.The statistics and the equilibration of the system are both better in this simulation set.

FIG. S17 .
FIG. S17.Parallel flow coupling (JC > 0): ground states and final configurations.(a) Example of a ground state configuration of our model when JC > 0. The two constraints are (i) the loops are fully-packed (blue-red-gray at each vertex) and (ii) all couplings are parallel (σiσj = +1 for each gray edge).They imply that the lowest energy configurations correspond to parallel streamlines that wind around the simulation box.(b) A configuration reached by our numerical simulation.It is not a ground state configuration of H.It however features a collection of nested loops.
FIG. S18.Parallel flow coupling (JC > 0): numerical results.(a) Gyration radius versus loop size.We uncover two regimes: for short loops (L < 100) the exponent is ν1 ≈ 0.85 while for large loops we find ν2 ≈ 0.60.(b) The correlation C(r) is consistent with a r −1 decay.(c) Distribution of loop sizes.