1. Introduction
Active turbulence, similar to its passive classical counterpart, is characterized by the emergence of highly complex bulk flow dynamics (Alert, Casademunt & Joanny Reference Alert, Casademunt and Joanny2022; Bentkamp et al. Reference Bentkamp, Drivas, Lalescu and Wilczek2022; Matsuzawa et al. Reference Matsuzawa, Mitchell, Perrard and Irvine2023). Active fluids based on motile bacteria (Sokolov et al. Reference Sokolov, Aranson, Kessler and Goldstein2007; Wensink et al. Reference Wensink, Dunkel, Heidenreich, Drescher, Goldstein, Löwen and Yeomans2012; Dunkel et al. Reference Dunkel, Heidenreich, Drescher, Wensink, Bär and Goldstein2013b), molecular motors (Sanchez et al. Reference Sanchez, Chen, DeCamp, Heymann and Dogic2012) and self-propelled colloids (Bricard et al. Reference Bricard, Caussin, Desreumaux, Dauchot and Bartolo2013) can display a rich set of topological structures, from spontaneously forming and annihilating point-defects in two-dimensional (2-D) films (Doostmohammadi et al. Reference Doostmohammadi, Adamer, Thampi and Yeomans2016) to entangled vortex lines in three-dimensional (3-D) bulk flows (Čopar et al. Reference Čopar, Aplinc, Kos, Žumer and Ravnik2019). Building on classic work on the statistical mechanics of point defects (Onsager Reference Onsager1949; Kosterlitz & Thouless Reference Kosterlitz and Thouless1973), the dynamics and statistics of topological defects have been extensively studied in (quasi) 2-D active fluids (Thampi, Golestanian & Yeomans Reference Thampi, Golestanian and Yeomans2014; Giomi Reference Giomi2015; James, Bos & Wilczek Reference James, Bos and Wilczek2018; Chardac et al. Reference Chardac, Hoffmann, Poupart, Giomi and Bartolo2021). By contrast, the diverse and complex singular structures realized by active flows in 3-D space (Binysh et al. Reference Binysh, Kos, Čopar, Ravnik and Alexander2020) were, until recently, inaccessible to experimental and numerical studies. With modern experimental imaging techniques (Duclos et al. Reference Duclos2020) and simulation methods, it is now possible to probe 3-D topological structures and their statistics (Kralj, Ravnik & Kos Reference Kralj, Ravnik and Kos2023).
Topological approaches have helped progress theoretical fluid mechanics since the observation by Moffatt (Reference Moffatt1969) that the tangling of vortex lines is related to the total helicity of the associated ideal incompressible flow. The total helicity 
$H$ is an inviscid invariant of incompressible flow defined by the integral (Moffatt Reference Moffatt1969)
\begin{equation} H(t) = \int_V \,\mathrm{d}^3\boldsymbol{x}\ \boldsymbol{u}(\boldsymbol{x},t)\boldsymbol{\cdot}\boldsymbol{\omega}(\boldsymbol{x},t), \end{equation}
where 
$\boldsymbol {u}$ is the flow velocity and 
$\boldsymbol {\omega } = \boldsymbol {\nabla } \times \boldsymbol {u}$ is its associated vorticity. Moffatt showed that for a system of 
$N$ closed and isolated vortex lines in simply connected domains, the total helicity can be expressed as
\begin{equation} H = \sum_{i,j=1}^N \kappa_i \kappa_j \mathcal{L}_{ij}, \end{equation}
where 
$\kappa _i$ is the circulation around the 
$i$th vortex tube, and 
$\mathcal {L}_{ij}$ is the linking number between the centrelines of the 
$i$th and 
$j$th vortex tubes, a topological measure counting the signed integer number of times the 
$j$th vortex tube wraps around the 
$i$th. This connection between helicity and flow topology has been measured in hydrodynamic experiments (Kleckner & Irvine Reference Kleckner and Irvine2013; Scheeler et al. Reference Scheeler, van Rees, Kedia, Kleckner and Irvine2017) and applies generally to solenoidal fields with discrete localization of their curl, including magnetic fields (Berger & Field Reference Berger and Field1984) or quantum flows (Hänninen & Baggaley Reference Hänninen and Baggaley2014; Zuccher & Ricca Reference Zuccher and Ricca2017).
Intuitively, 
$H$ measures the winding of vortex lines around each other, and non-zero helicities indicate chiral flows in which vortex lines curl around each other in a preferred orientation. While this connection between helicity and topology is far-reaching, in most flows of interest, the topology of the vorticity field is much more complex than a set of potentially interlinked discrete vortex tubes. As turbulence emerges, experimental resolution of vortex lines becomes more challenging in the absence of localized vortex tubes (Matsuzawa et al. Reference Matsuzawa, Mitchell, Perrard and Irvine2023), and (1.2) does not directly hold: for general vorticity fields, vortex lines are not necessarily closed, and an asymptotic formulation of (1.2) due to Arnold (Reference Arnold1974) provides an extension of Moffatt's result. Although theoretically useful, Arnold's formulation has been relatively underutilized in practice due to the high computational cost involved in the computation of the relevant topological quantities. It is thus desirable to establish practical means to characterize the topological structure of diffuse vorticity fields, independently of the dynamics through which these arise.
Here, building on these core ideas, we study the statistics of vortex line entanglement in a turbulent flow governed by the incompressible generalized Navier–Stokes (GNS) equations (Słomka & Dunkel Reference Słomka and Dunkel2017b; Słomka, Suwara & Dunkel Reference Słomka, Suwara and Dunkel2018; Supekar et al. Reference Supekar, Heinonen, Burns and Dunkel2020)
$$\begin{gather} \partial_t \boldsymbol{u} + (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} =-\boldsymbol{\nabla} p + \varGamma_0 \nabla^2 \boldsymbol{u}-\varGamma_2 \nabla^4 \boldsymbol{u} + \varGamma_4 \nabla^6 \boldsymbol{u}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$which model flows with advected active constituents driving a generic linear instability (Rothman Reference Rothman1989; Beresnev & Nikolaevskiy Reference Beresnev and Nikolaevskiy1993; Tribelsky & Tsuboi Reference Tribelsky and Tsuboi1996). The GNS equations (1.3) provide a coarse-grained phenomenological description of active fluid flows on scales larger than those of the individual microorganisms or self-propelled particles that drive these flows. Linkmann et al. (Reference Linkmann, Marchetti, Boffetta and Eckhardt2020) showed how such GNS equations can be derived from a more complex active suspension model.
Nonlinear advection is nominally negligible in the low-Reynolds-number regime typical of microfluidic experiments with passive fluids or dilute active suspensions. By contrast, the presence of active stresses in non-dilute microswimmer suspensions can significantly alter the viscous balance, effectively cancelling the fluid viscosity (López et al. Reference López, Gachelin, Douarche, Auradou and Clément2015; Nambiar, Nott & Subramanian Reference Nambiar, Nott and Subramanian2017) and creating an unstable band of modes. These modes can saturate nonlinearly via the advective term, and the resulting system can therefore have a large effective Reynolds number and exhibit turbulent-like behaviour (Dunkel et al. Reference Dunkel, Heidenreich, Bär and Goldstein2013a; Koch & Wilczek Reference Koch and Wilczek2021). The GNS equations model this behaviour by inducing a band-limited linear instability via the 
$\varGamma$ terms which saturates in a finite-amplitude, statistically stationary flow. For 
$\varGamma _0, \varGamma _4 >0$ and 
$\varGamma _2 < 0$, the parameters 
$(\varGamma _0, \varGamma _2, \varGamma _4)$ together define a characteristic energy injection length scale 
$\varLambda$ and bandwidth 
$\kappa$, along with a characteristic time scale 
$\tau\kern-0.8pt$:
\begin{equation} \varLambda = {\rm \pi}\sqrt{-\frac{2\varGamma_4}{\varGamma_2}}, \quad \kappa = \left(-\frac{\varGamma_2}{\varGamma_4} - 2\sqrt{\frac{\varGamma_0}{\varGamma_4}}\right)^{1/2}, \quad \tau = \frac{2\varGamma_4}{\varGamma_2} \left(\varGamma_0 - \frac{\varGamma_2^2}{4\varGamma_4} \right)^{-1}. \end{equation}
Physically, 
$\varLambda$ prescribes the characteristic vortex size, while 
$\kappa$ determines how many modes with wavelengths around 
$\varLambda$ are excited by the linear instability. The phenomenology of typical microfluidic experiments with bacterial suspensions (Dunkel et al. Reference Dunkel, Heidenreich, Drescher, Wensink, Bär and Goldstein2013b; Wioland et al. Reference Wioland, Woodhouse, Dunkel and Goldstein2016; Čopar et al. Reference Čopar, Aplinc, Kos, Žumer and Ravnik2019; Peng, Liu & Cheng Reference Peng, Liu and Cheng2021) can be recovered by setting 
$\varLambda \sim 75$ 
${\rm \mu}$m, 
$\kappa \sim 30$ mm
$^{-1}$ and with a characteristic speed 
$U = 2{\rm \pi} \varLambda /\tau \sim 10\,\mathrm {\mu }{\rm m}\ {\rm s}^{-1}$ (Słomka & Dunkel Reference Słomka and Dunkel2017b). Other, more complex linear dispersion relations produce structurally similar flow patterns (Słomka et al. Reference Słomka, Suwara and Dunkel2018).
In this active fluid context, the dynamics of helicity become especially interesting: the combination of a linear instability and nonlinear advection can lead to a partial 3-D inverse energy transport mediated by chiral symmetry breaking (Biferale, Musacchio & Toschi Reference Biferale, Musacchio and Toschi2012; Słomka & Dunkel Reference Słomka and Dunkel2017b). Following the traditional characterization of the inverse energy transport in 2-D passive turbulence through vortex coarsening (Boffetta & Ecke Reference Boffetta and Ecke2012), here we investigate the 3-D topology of vortex lines in 3-D GNS flows.
While previous analytical and numerical studies of (1.3) and other active fluids in 3-D have focused on unbounded domains (Słomka & Dunkel Reference Słomka and Dunkel2017b; Urzay, Doostmohammadi & Yeomans Reference Urzay, Doostmohammadi and Yeomans2017; Słomka et al. Reference Słomka, Suwara and Dunkel2018), analysing topological structures in a typical experimental set-up requires confined simulations. Moreover, Moffatt's and Arnold's theorems, like many theoretical tools relying on topological properties of the ambient space, are technically applicable only in simply connected domains. We hence choose a 3-D ball of radius 
$R$ as our simulation domain, which is qualitatively similar to experimentally realisable microfluid cavities (Wioland et al. Reference Wioland, Woodhouse, Dunkel, Kessler and Goldstein2013).
In this work, we leverage spectral direct numerical simulations (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020) and recent methods to efficiently compute linking numbers (Qu & James Reference Qu and James2021) to explore the topological structure of a numerical realization of active turbulence in confinement (figure 1). However, the methodology developed here is applicable to any incompressible flow, including passive and active fluids. Our active model spontaneously breaks spatial parity and, consistent with earlier results for periodic domains (Słomka & Dunkel Reference Słomka and Dunkel2017b), produces a quasi-Beltrami flow in the bulk of a closed domain, as will be shown in § 2. We then characterize the topological structure of the emergent Beltrami flow by numerically computing the entanglement statistics of vortex lines in § 3. To validate this characterization, we observe that the mean linking number between two vortex lines converges to the asymptotic results due to Arnold (Reference Arnold1974). Beyond this result, the full distribution of linking numbers is well described by a k-Gamma distribution, in agreement with an entropic argument previously encountered in granular and living matter for distributions of one-sided random variables with constraints on the mean (Aste & Di Matteo Reference Aste and Di Matteo2008; Atia et al. Reference Atia2018; Day et al. Reference Day2022). This statistical argument is detailed in § 4.
Figure 1. Simulations of the active GNS model (1.3) confined in the ball demonstrate spontaneous parity-breaking and helical flows. Structures are more pronounced for a narrow spectral energy injection bandwidth ((a) 
$\kappa =4/R$) than a wider one ((b) 
$\kappa =16/R$). Iso-value surfaces of the helicity density 
$h=\boldsymbol {u}\boldsymbol{\cdot} \boldsymbol {\omega }$ are shown, rescaled by the standard deviation 
$\sigma (h)$. Regions of positive helicity (red) dominate in the bulk and negative regions (blue) are mostly present near the boundary. The vortex scale parameter is 
$\varLambda = R/8$ in both simulations, and time units are chosen such that 
$\tau =1$.
2. Chiral symmetry breaking in linearly forced active flows
The unbounded GNS equations admit exact chiral solutions. These solutions are Beltrami flows, in which the vorticity is colinear with the velocity (Słomka & Dunkel Reference Słomka and Dunkel2017b). More specifically, the exact GNS solutions have a velocity field which is an eigenfunction of the curl operator
\begin{equation} \boldsymbol{\omega} =\boldsymbol{\nabla} \times \boldsymbol{u} = \lambda \boldsymbol{u}, \end{equation}
with 
$\lambda$ the eigenvalue corresponding to the characteristic wavenumber of the mode; such solutions are sometimes called Trkalian flows (Lakhtakia Reference Lakhtakia1994). The initial linear instability drives modes with 
$\lambda \approx {\rm \pi}/\varLambda$ and both positive and negative helicities. Triadic interactions break the symmetry and spontaneously select an overall handedness for the flow (Słomka et al. Reference Słomka, Suwara and Dunkel2018). Simulations in periodic domains that start from initially random small velocity fields therefore spontaneously produce statistically stationary chiral flows with non-zero mean helicity. However, it has remained an open question whether such solutions can robustly manifest themselves in the presence of boundaries.
Using the spectral solver Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), we simulate a confined GNS flow inside the three-dimensional ball of radius 
$R$ (figure 1). Starting from an initially small random flow, we evolve the GNS equations (1.3) subject to the boundary conditions:
$$\begin{gather} \boldsymbol{u}(\boldsymbol{R},t) = 0, \end{gather}$$
$$\begin{gather}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}(\boldsymbol{R},t) = 0, \end{gather}$$
$$\begin{gather}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\nabla}\boldsymbol{u}(\boldsymbol{R},t) = 0, \end{gather}$$
where 
$\boldsymbol {n}$ is the outward pointing unit normal vector and 
$\boldsymbol {R}$ any point on the boundary. The no-slip condition here ensures there is no normal vorticity at the boundary (a necessary condition for Arnold's theorem), and the higher order terms are chosen to simply suppress shear at the boundary, but other choices are possible (Słomka & Dunkel Reference Słomka and Dunkel2017a). Throughout, we set 
$\varLambda = R/8$ and adopt time units such that 
$\tau =1$. To simulate the GNS equations (1.3) with boundary conditions (2.2), we discretize the ball along the 
$(r,\theta, \phi )$ coordinates using 
$(N_r, N_\theta, N_\phi ) = (64, 64, 128)$ grid points. Time stepping is done using a third-order four-stage implicit–explicit Runge–Kutta scheme (Ascher, Ruuth & Spiteri Reference Ascher, Ruuth and Spiteri1997). Additional details of the numerical implementation and initial condition construction are summarized in the Appendix.
As has been observed in periodic domains (Słomka & Dunkel Reference Słomka and Dunkel2017b), after an initial transient, the GNS dynamics lead to the spontaneous emergence of vortical flows which saturate at a finite energy 
$E(t) = (1/2)\int _V \,\mathrm {d}^3\boldsymbol {x} |\boldsymbol {u}(\boldsymbol {x},t)|^2$ (figures 1 and 2). The asymmetry between positive and negative helicity density regions suggests that the emergent flow violates parity invariance (figure 1): while the GNS equations – including the chosen boundary conditions – are invariant under the transformation 
$\boldsymbol {x} \rightarrow -\boldsymbol {x}$, solutions with a non-zero total helicity are not invariant under this transformation. To quantify the extent of this parity-symmetry breaking, we compute the helicity of the flow 
$H(t)$ through its integral definition. We find that 
$H(t)$, like the energy 
$E(t)$, starts at a small value and grows until saturating at a much larger amplitude. However, unlike the energy, the steady-state helicity has a sign that is determined by the random initial condition (figure 2). It is interesting to note the variability in transient behaviour between samples, with the presence of dynamics on multiple time scales. This phenomenon is reminiscent of mode competition in nonlinear systems such as multimode lasers (Hodges et al. Reference Hodges, Munroe, Gadomski, Cooper and Raymer1997) and population dynamics (Hastings et al. Reference Hastings, Abbott, Cuddington, Francis, Gellner, Lai, Morozov, Petrovskii, Scranton and Zeeman2018; Morozov et al. Reference Morozov, Abbott, Cuddington, Francis, Gellner, Hastings, Lai, Petrovskii, Scranton and Zeeman2020).
Figure 2. Spontaneous chiral flows in confined GNS. Energy and total helicity as a function of time for different initial conditions: (a) 20 simulations with 
$\varLambda = R/8, \kappa =4/R$; (b) 19 simulations with 
$\varLambda = R/8, \kappa =16/R$. In both panels (a,b), the initial conditions spontaneously break parity symmetry, and 
$\tau =1$ sets the time unit.
What is the structure of these emergent chiral solutions? To compare our flows with the expected Beltrami solutions in the bulk, we compute a ‘Beltrami factor’ 
$\beta = \boldsymbol {u}\boldsymbol{\cdot} \boldsymbol {\omega }/ (\lambda |\boldsymbol {u}|^2)$ with 
$\lambda = {\rm \pi}/\varLambda$. If the flow followed the structure of the periodic solutions (2.1), we would expect this measure to be peaked around 
$\pm 1$ in the bulk, with 
$+1$ corresponding to positive helicity solutions. Indeed, as we decrease 
$\kappa$ and fewer modes are excited, 
$\beta$ peaks around 
$1$, with the notable appearance of secondary peaks (figure 3a). Those new peaks can be simply explained: in confined domains, solutions to (2.1) cannot also be solutions to the GNS equations (1.3) subject to the chosen boundary conditions (2.2). This frustration leads to the appearance of a boundary layer (as can be seen in figure 1 for 
$\kappa =4/R$), and indeed 
$\beta$ is peaked around unity in the bulk of the sphere (figure 3b). By Taylor expanding near 
$r=R$ to third order with the chosen boundary conditions (2.2), and matching to the typical third derivative in the bulk 
${\sim }U/\varLambda ^3$, the characteristic boundary layer scale is expected to be 
$w = 6^{1/3}\varLambda \approx 1.8 \varLambda$, which matches well with our simulations (figure 3b,c).
Figure 3. GNS spontaneously produces bulk quasi-Beltrami flows for narrow unstable bandwidth. (a) Beltrami factor 
$\beta = \boldsymbol {u}\boldsymbol{\cdot} \boldsymbol {\omega }/ (\lambda |\boldsymbol {u}^2|)$ with 
$\lambda = {\rm \pi}/\varLambda$ is increasingly peaked for systems with narrower spectral bandwidth 
$\kappa$. The depicted flows all have positive helicity, but opposite parity solutions appear with equal probability. (b) Two-dimensional histogram of the Beltrami factor against the radial position, revealing an approximate Beltrami flow in the bulk of the ball with adjustments in a boundary layer of relative thickness 
$\hat {w} = 6^{1/3}\varLambda /R$. At higher bandwidths (inset), more modes are excited and the Beltrami factor is less clustered. In all simulations in panels (a,b), 
$\varLambda = R/8$. (c) Beltrami factor for simulations with 
$\kappa R=4$ and varying 
$\varLambda$ support the expected boundary layer scaling. The histograms in panels (b,c) are constructed from 
$10^4$ uniformly random sample points. All values are in simulation units where 
$R=1, \tau = 1$.
The generalized Navier–Stokes equation hence spontaneously generate quasi-Beltrami flows in the bulk for narrow energy injection bandwidths. Why do GNS solutions converge to such flows? Chiral symmetry breaking has been explained in previous work by noticing that the advection term selects for chiral solutions in the bulk in the presence of energy injection by linear instability (Słomka et al. Reference Słomka, Suwara and Dunkel2018). This selection effect is theorized to be more pronounced as the energy injection bandwidth 
$\kappa$ narrows, in accordance with our numerical observations, as both the absolute value of the helicity and 
$\beta$ decrease with larger bandwidths (figures 2b and 3). As Beltrami flows minimize enstrophy 
$\mathcal {E} = (1/2)\int \,\mathrm {d}^3\boldsymbol {x} \boldsymbol {\omega }^2$ at fixed helicity (Woltjer Reference Woltjer1958) and are steady solutions to the Euler equations, one might speculate that in our effectively low-viscosity flows (López et al. Reference López, Gachelin, Douarche, Auradou and Clément2015), the selection of helical modes naturally leads to such Beltrami flows away from boundaries (Słomka & Dunkel Reference Słomka and Dunkel2017b). While out-of-equilibrium dynamics do not necessarily follow any extremization principle, we note that this property of Beltrami flows is purely geometric and does not depend on the nature of the flow. It would thus be interesting to see under what conditions other helical turbulence generation mechanisms also produce Beltrami flows, and whether this regime could be realized in microfluidic experiments using semi-dense bacterial suspensions (Wioland et al. Reference Wioland, Woodhouse, Dunkel, Kessler and Goldstein2013) or other biological or synthetic active matter.
3. Quantifying chiral symmetry breaking through vortex linking statistics
In the previous section, we showed that the GNS flow in the ball spontaneously produces quasi-Beltrami flows with non-zero helicity. To connect the helicity of the flow to the linking statistics of vortex lines, we follow an approach inspired by theorems by Moffatt (Reference Moffatt1969) and Arnold (Reference Arnold1974).
Both theorems are concerned with the instantaneous entanglement of vortex lines, which are defined as the streamlines (in the mathematical sense) of the vorticity field, satisfying the differential equation
\begin{equation} \frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}\tilde{s}} = \boldsymbol{\omega}(\boldsymbol{r}(\tilde{s})) , \end{equation}
where 
$\tilde {s}$ is a parameter with units of length 
$\times$ time. Here and in what follows, this integration is always performed at a fixed simulation time point 
$t$, considering the vorticity field 
$\boldsymbol {\omega }$ as a ‘frozen-in’ structure at time 
$t$. To avoid numerical issues due to varying magnitude of 
$\boldsymbol {\omega }$, we consider the equivalent differential equation
\begin{equation} \frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}s} = \frac{ \boldsymbol{\omega}}{| \boldsymbol{\omega}|}(\boldsymbol{r}(s)) \end{equation}
now re-parametrized such that a vortex line integrated from an initial position 
$\boldsymbol {r}(0)=\boldsymbol {r}_0$ over 
$s\in [0, L]$ has length 
$L$. We numerically integrate (3.2) by using a linear interpolation of the field 
$\boldsymbol {\omega }$ between the quadrature nodes of the spectral direct numerical simulation.
Vortex lines are experimentally accessible flow structures in the limit of localized vorticity, as buoyant particles such as bubbles in water are attracted to regions of high vorticity while heavier particles are expelled away from them. This makes vortex lines readily visible in flows where vorticity is confined to tube-like regions (Durham et al. Reference Durham, Climent, Barry, De Lillo, Boffetta, Cencini and Stocker2013; Kleckner & Irvine Reference Kleckner and Irvine2013). However, vortex line geometry is often complex even in the presence of relatively simple flow fields, as illustrated by the chaotic field lines in the Arnold–Beltrami–Childress (ABC) flow (Dombre et al. Reference Dombre, Frisch, Greene, Hénon, Mehr and Soward1986; Qin & Liao Reference Qin and Liao2023). Integration of vortex lines in our helical flows indeed lead to erratic trajectories, where initially close vortex lines rapidly diverge and tangle around each other (figure 4).
Figure 4. Vortex lines, corresponding to integral curves of the vorticity field at fixed time 
$t$ as defined in (3.1), starting at evenly spaced intervals along the black line, rapidly diverge and tangle around each other. (a) Lines of length 
$L\approx 0.3R$; the two longer lines have 
$L\approx 0.9R$. Colour indicates starting position. (b) The same lines as in panel (a) extended to 
$L=5R$. The vorticity field is from the statistically stationary turbulent state with 
$\varLambda = R/8$, 
$\kappa =4/R$.
To measure pair-wise entanglement of vortex lines, we use the linking number between two oriented closed curves 
$\gamma _1$ and 
$\gamma _2$ defined by Gauss’ integral formula (Qu & James Reference Qu and James2021)
\begin{equation} \mathcal{L}(\gamma_1, \gamma_2) = \frac{1}{4{\rm \pi}}\oint_{\gamma_1} \oint_{\gamma_2} \frac{(\boldsymbol{r}_1 - \boldsymbol{r}_2)}{|\boldsymbol{r}_1 - \boldsymbol{r}_2|^3}\boldsymbol{\cdot} (\mathrm{d}\boldsymbol{r}_1 \times \mathrm{d}\boldsymbol{r}_2).\end{equation}
This integer-valued quantity counts the signed number of times one curves winds around the other. The linking number possesses notable properties: if one curve is reversed, the linking number flips sign and it is symmetric since 
$\mathcal {L}(\gamma _1, \gamma _2) = \mathcal {L}(\gamma _2, \gamma _1)$. Most importantly, 
$\mathcal {L}$ is invariant under continuous deformation of the curves. The linking number is therefore a topological invariant playing a central role in the study of knots and linked curves (Vologodskii & Cozzarelli Reference Vologodskii and Cozzarelli1994; Kauffman Reference Kauffman1995; Panagiotou Reference Panagiotou2019). An equivalent viewpoint is to define 
$\mathcal {L}$ as the sum of half-integer contributions from each crossing of the curves under a planar projection, with the sign depending on the relative orientation of the curves (figure 5a). As an invariant, 
$\mathcal {L}$ is independent of the choice of the projection plane.
Figure 5. Linking numbers are topological invariants measuring the entanglement of oriented closed curves. (a) Example configurations with their linking numbers. Positive and negative crossings are illustrated by 
$\oplus$ and 
$\ominus$. (b) Linking number of two example open curves closed using the scheme used to apply Arnold's theorem. (c) Illustration of the additive property of linking numbers with respect to curve concatenation. (d) The additive property can be used to compute linking numbers as sums of contributions from subsegments of two curves.
While the definition of 
$\mathcal {L}$ in (3.3) calls for closed curves, vortex lines in complex flows are not closed in general. To leverage the connections between linking numbers and helicity, we thus have to consider slightly modified vortex lines. Consider two open vortex lines 
$\gamma _1$ and 
$\gamma _2$ of length 
$L$ starting at points 
$\boldsymbol {x}_1$ and 
$\boldsymbol {x}_2$, respectively. We then define the asymptotic linking number 
$\varLambda _L(\boldsymbol {x}_1,\boldsymbol {x}_2)$ between 
$\gamma _1$ and 
$\gamma _2$ as follows. We construct the curve 
$\tilde {\gamma }_i$ as by closing the curve 
$\gamma _i$ by a straight segment connecting its start and end points. We then note 
$\varLambda _L(\boldsymbol {x}_1,\boldsymbol {x}_2)$ the linking number of the two curves 
$\tilde {\gamma }_1$, 
$\tilde {\gamma }_2$ (figure 5b) normalized by 
$T_i$ and 
$T_j$:
\begin{equation} \varLambda_L(\boldsymbol{x}_1,\boldsymbol{x}_2) = \frac{1}{T_1 T_2}\mathcal{L}(\tilde{\gamma}_1, \tilde{\gamma}_2), \end{equation}
where 
$T_i = \int _{\gamma _i} \,\mathrm {d}s/|\omega (\boldsymbol {r}(s))|$ is the integrated inverse circulation, which has units of length 
$\times$ time; the product 
$T_iT_j$ is the normalization factor of the asymptotic linking number.
With these preliminary definitions, Arnold (Reference Arnold1974) provides a connection between the helicity of an incompressible flow and the asymptotic linking number of vortex lines via the asymptotic equality:
\begin{equation} H=\int_V \,\mathrm{d}^3\boldsymbol{x}\ (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\omega})(\boldsymbol{x}) = \frac{1}{V^2}\iint_{V\times V}\,\mathrm{d}^3\boldsymbol{x}_1\,\mathrm{d}^3\boldsymbol{x}_2 \lim_{L\rightarrow+\infty}\varLambda_L(\boldsymbol{x}_1, \boldsymbol{x}_2). \end{equation}
This equality between the helicity and the volume averages of the asymptotic linking numbers is valid on any simply connected domain for any incompressible velocity field, as long as the normal component of the vorticity vanishes at the boundary (
$\boldsymbol {n}\boldsymbol{\cdot} [\boldsymbol {\nabla } \times \boldsymbol {u}] = 0$). Since our simulations have a no-slip boundary condition 
$\boldsymbol {u}=0$, by Stokes’ theorem, the vorticity flux through any arbitrary loop drawn on the boundary cancels, and we can apply (3.5).
To characterize the topology of the active Beltrami flow, we hence construct a statistical ensemble of linking numbers using the numerical solutions of (3.2), verifying our approach against the expectation of (3.5) for the mean of this distribution.
As a first step towards an estimate of the helicity using Arnold's theorem, we integrate 
$N$ vortex lines of length 
$L$ and then close them as in figure 5(b). To compute the linking number between two such closed vortex lines, which are numerically represented as a set of connected line segments with lengths 
$\delta <10^{-2}R$, a naive discretization of the Gauss integral of (3.3) between two lines of length 
$L$ would require 
$O(L^2/\delta ^2)$ operations. This quadratic scaling leads to a prohibitive computational cost for large linking numbers. However, 
$N$-body simulations often require the evaluation of integrands which decay as 
$\propto 1/r^2$, as is the case for (3.3). For this class of functions, one can leverage Barnes–Hut and fast multipole-type methods to bring down the algorithmic complexity of the linking number computation to 
$O((L/\delta )\log (L/\delta ))$. Qu & James (Reference Qu and James2021) designed and implemented such methods for computing linking numbers, along with topology-preserving curve simplification algorithms in a publicly available C++ package. We have built and released a Python wrapper for their code (see Data availability statement).
To monitor convergence of the mean of the asymptotic linking numbers to 
$H$ as a function of line length 
$L$, we can exploit the linearity of the curve integrals in (3.3). By decomposing the linking number into contribution from subloops, we can avoid recomputing linking numbers from scratch (Moffatt Reference Moffatt1969). Let 
$\gamma _1 \gamma _2$ denote the concatenation of the oriented closed curves 
$\gamma _1$ and 
$\gamma _2$ sharing a start and end point, then we have 
$\mathcal {L}(\gamma _1\gamma _2,\gamma _3) = \mathcal {L}(\gamma _1,\gamma _3) + \mathcal {L}(\gamma _2,\gamma _3)$ (figure 5c). Contributions from subloops can then be summed up to recover the full linking number of longer curves (figure 5d).
With these ingredients in hand, we can compute the linking number distribution and construct an ‘Arnold estimate’ of the helicity, by constructing a Monte Carlo approximation to the integral in the right-hand side of (3.5) as follows.
(i) Sample
$N$ initial points $\{\boldsymbol {x}_i\}$ uniformly in the domain.(ii) Integrate vortex lines and their inverse circulation for a length
$L$. Closing the vortex lines by a line segment, we obtain a set of $\{g(\boldsymbol {x}_i), T_i\}$.(iii) Compute the
$N(N-1)/2$ distinct linking numbers $\lambda _{ij} := \mathcal {L}(g(\boldsymbol {x}_i),g(\boldsymbol {x}_j))$ using fast methods (Qu & James Reference Qu and James2021).(iv) Normalize the linking numbers to obtain the approximations
$\varLambda _L$ and average those contributions to estimate the average helicity :
(3.6)\begin{equation} \hat{H} = \frac{2}{N(N-1)}\sum_{0< i< j\leq N} \frac{1}{T_i T_j}\lambda_{ij}. \end{equation}
Applying this program, we study the convergence of the estimate 
$\hat {H}$ to 
$H$ as a function of 
$N$ and 
$L$ for a snapshot of a given simulation at a fixed time point once the GNS has entered the statistically stationary state.
Integrating vortex lines, we find normalization factors 
$T_iT_j$ increasing as 
${\sim } L^2$, as can be expected for long vortex lines with 
$L\gg R$ which traverse the entire domain such that 
$T_i \approx L/\langle |\boldsymbol {\omega }|\rangle$, with 
$\langle. \rangle$ denoting the volume average (figure 6a). Note that boundaries are irrelevant to the computation: as 
$T_i$ diverges to 
$+\infty$ when 
$\boldsymbol {u}$ vanishes, vortex lines intersecting the boundary layer have very large normalization factors. The contributions from such ‘boundary’ vortex lines are hence suppressed from the estimate (3.6), and the mean of 
$T_iT_j$ is larger than its mode. Once the vortex lines are computed, we compute our ensemble of linking numbers; consistent with the implication of (3.6) that 
$\varLambda _L$ must converge to a finite value, our computed linking numbers have their average value scaling with 
$L^2$ (figure 6b). It is interesting to note that vortex lines longer than a few domain radius will almost certainly link with other vortex lines with the sign of the total helicity. Additionally, the linking number between vortex lines is independent of vortex line starting position; this is consistent with the picture that vortex line integration is chaotic (figure 6b, inset).
Figure 6. The average linking number converges to the total helicity as vortex line length is increased. (a) Heatmap of the normalization factor 
$T_iT_j$ as a function of vortex line length 
$L$ averaged over 
$N = 1000$ vortex lines, in (a i) linear and (a ii) log scale. Continuous lines indicate mean, dashed lines standard deviation. Notice that the mean sits above the mode of the distribution, as the normalization factor distribution has a long tail due to interactions with the boundaries. (b) Heatmap of the pair-wise linking number 
$\mathcal {L}$ as a function of vortex line length 
$L$, in (b i) linear and (b ii) log scale for a positive helicity flow. For better visualization, each column in panels (a,b) is normalized to the column-wise standard deviation. Results are averaged over 
$N = 1000$ vortex lines and their 
$N(N-1)/2 = 5\times 10^5$ pairs. The logarithmic plot shows an asymptotic scaling of 
$\mathcal {L} \sim L^2$. Inset shows the set of linking numbers 
$\mathcal {L}$ for 
$L=100R$ as a function of the distance 
$d$ between the vortex lines starting points. There is no correlation between 
$\mathcal {L}$ and 
$d$. (c) Helicity estimate as a function of vortex line length for increasing number of vortex lines. (d) Helicity estimate for various samples (different colours) computed from 
$N=500$ vortex lines as a function of 
$L$. Remarkably, even short lines capture the helicity sign, suggesting that even limited observations of vortex lines could be used to detect chiral symmetry breaking in experiments.
Combining linking numbers and normalization factors to compute 
$\hat {H}$, we find rapid convergence with length and number: 
$N=125$ vortex lines of length 
${\sim }30R$ give an estimate accurate within 10 %. Remarkably, even short vortex lines lead to the correct sign and order of magnitude of the helicity, across tested samples (figure 6d).
Settling on 
$N=150, L=100R$ to construct our Arnold estimate 
$\hat {H}$, we run the above algorithm at various time points to monitor the time-evolution of vorticity linking (figure 7). The Arnold estimate is an accurate estimate of the helicity during the initial period, linear instability and saturation phase, with an approximately constant relative error (figure 7a,b). In line with our previous observation that even short vortex lines led to the correct helicity sign, 
$\hat {H}$ has the correct sign even when the helicity magnitude is close to zero. This raises the interesting possibility of tracking chiral symmetry breaking in experiments using tracer particles to uncover vortex lines, especially at low tracer densities that would make standard velocity reconstruction methods challenging.
Figure 7. Tracking chiral symmetry-breaking through vorticity linking. (a) Direct integration of the helicity density (solid line) matches the Arnold estimate of the helicity by vortex linking number statistics (solid circles). Dashed line indicates the average helicity at steady state. (b) The Arnold estimate is accurate even at short times, and notably captures the helicity sign change, with red dots denoting negative values and blue dots indicating positive values. The 
$y$-axis is linear between 
$[-10^{-7}, 10^{-7}]$ and logarithmic elsewhere. (c) Time evolution of the probability distribution of pair-wise linking numbers. Colour indicates time. (d) Mean, standard deviation and moment coefficient of skewness of the linking number distribution as a function of time. The final non-zero mean reflects chirality of flow, while increasing skewness indicates the departure from Gaussian statistics. We note that at time 
$10$, an outlier point with skewness 
$40$ is not shown. Results are shown for 
$N= 150$ vortex lines of length 
$L=100R$. Colour of scatter points in panels (a,c,d) all follow the colour map of panel (c).
The dynamics of vortex lines is key to understanding the emergence of fine structures in turbulence through their connection to helicity as an inviscid invariant (Scheeler et al. Reference Scheeler, van Rees, Kedia, Kleckner and Irvine2017; McKeown et al. Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2020; Matsuzawa et al. Reference Matsuzawa, Mitchell, Perrard and Irvine2023). In active flows, however, helicity can be created and destroyed. Our techniques allow the characterization of the statistically averaged topology of the flow as it evolves to a statistically stationary state with non-zero total helicity. Here, our construction of 
$\hat {H}$ provides us with the full distribution 
$p(\mathcal {L})$ of linking numbers as a function of time, beyond Arnold's result on the mean degree of linkage of vortex lines (figure 7c). Studying the moments of this distribution, one finds that, as expected, the mean linking number goes from 
$0$ to finite non-zero number reflecting the emergence of chiral flows (figure 7d). The behaviour of the standard deviation and skewness are however non-trivial: the standard deviation shows large fluctuations during the instability growth phase, while the distribution displays non-zero skewness at late times. In the next section, we will use general constraints on the linking number distribution to rationalize its statistics at steady state.
4. The distribution of linking numbers obeys a maximum entropy law
The construction from the previous part allows us to obtain the distribution of pair-wise linking numbers of 
$N$ vortex lines of length 
$L$, which contains several notable features at steady state. First, for the strongly chiral flows considered here, all pairs of sufficiently long vortex lines are linked with probability 
$1$. Second, the mean 
$\langle \mathcal {L}\rangle$ is approximately constrained by Arnold's theorem to be equal to the flow's total helicity. Third, we find no correlation between linking numbers and vortex line starting points for sufficiently long vortex lines, suggesting ‘chaotic tangling’ and a notion of ergodicity in the system, with two randomly selected vortex lines eventually capturing the global helicity as their lengths tend to infinity. Together, these features suggest that in this geometrically and topologically complex system, statistical principles could explain the observed linking distribution.
Maximal-entropy reasoning has been successfully applied to explain the packing statistics of confined granular and living matter (Edwards & Oakeshott Reference Edwards and Oakeshott1989; Bi et al. Reference Bi, Henkes, Daniels and Chakraborty2015; Atia et al. Reference Atia2018; Day et al. Reference Day2022) and topological defect distributions in two-dimensional turbulence (Eyink & Sreenivasan Reference Eyink and Sreenivasan2006; Giomi Reference Giomi2015). As these problems naturally share features with our system of confined topological defects, the maximal-entropy method is a viable candidate to explain linking number statistics.
To apply a maximum entropy approach, we translate the above observations into constraints that the distribution of linking number must plausibly satisfy. The first observation implies that for a given length 
$L$, the linking numbers must be bounded from below 
$\mathcal {L} \geq \mathcal {L}_{min}$ for a positive helicity flow; in negative helicity flows, 
$\mathcal {L} \leq \mathcal {L}_{max}$. The second observation implies that the sum of linking numbers must be approximately equal to the flow helicity by (3.5):
\begin{equation} \sum_{i,j} \mathcal{L}_{ij} \approx \left(\frac{VNL}{\langle | \boldsymbol{\omega}|\rangle}\right)^2 H \equiv \bar{H}. \end{equation}
Here, we assumed long enough vortex lines 
$L \gg R$ such that we can take the flow to be homogeneous and 
$T_i\approx \langle | \boldsymbol {\omega }|\rangle /L$. Numerically, we do observe 
$\sum _{i,j} \mathcal {L}_{ij}$ to scale with 
$L^2$ (figure 6b). Finally, the third observation above suggests that there is effectively no correlation between linking numbers, or even perhaps between linking numbers of ‘long enough’ vortex line sub-segments; under this assumption, one can consider the linking number distribution as drawn from an emergent thermodynamic ensemble.
To proceed, we consider the distribution 
$p(\mathcal {L})$ which maximizes the Shannon entropy subject to the constraints outlined above. To this end, we consider 
$p(\mathcal {L})$ as the probability of the ‘macroscopic state’ where one vortex line of length 
$L$ links 
$\mathcal {L}$ times with another vortex line. Many ‘microscopic states’ corresponding to possible vortex line conformations are compatible with such macro-states. Following Aste & Di Matteo (Reference Aste and Di Matteo2008), we consider the Shannon entropy 
$\mathcal {S}$ written as
\begin{equation} \mathcal{S} = \sum_{\mathcal{L} \in \mathbb{Z}} - p(\mathcal{L})\log p(\mathcal{L}) + \sum_{\mathcal{L} \in \mathbb{Z}} p(\mathcal{L}) S(\mathcal{L}), \end{equation}
where 
$S(\mathcal {L})$ is the entropy of the state with linking 
$\mathcal {L}$. Under the assumption that all microscopic states are equiprobable, 
$S(\mathcal {L}) = \log \varOmega (\mathcal {L})$ with 
$\varOmega (\mathcal {L})$ the number of micro-states with linking 
$\mathcal {L}$. Under the maximum entropy principle, 
$p(\mathcal {L})$ is given by optimizing the entropy functional 
$\mathcal {S}$ under the helicity constraint
\begin{equation} \frac{1}{N^2}\bar{H} = \sum_{\mathcal{L} \in \mathbb{Z}} p(\mathcal{L}) \mathcal{L}. \end{equation}The solution of this optimization problem is given by a Boltzmann-type distribution:
\begin{equation} p(\mathcal{L}) \propto \varOmega(\mathcal{L}) {\rm e}^{-\mathcal{L}/\chi} \end{equation}
with 
$\chi ^{-1}$ the Lagrange multiplier fixing the helicity constraint. To fully determine the maximal-entropy distribution, the last step is to compute 
$\varOmega (\mathcal {L})$.
Motivated by our observation that even short vortex lines almost certainly link with others, we consider dividing a vortex line into 
$k$ sub-loops of an approximately constant size 
$\delta L$ which we see as characteristic of the domain, such that 
$k\approx L/\delta L$. A mesoscopic description of the linking of two vortex lines can be given by the linking numbers of each sub-loop of the first vortex line with the entire other line 
$\{\ell _i\}_{i=1,\ldots,k}$; in the case of a positive helicity flow, each sub-loop must link at least 
$\ell _{min} = \mathcal {L}_{min}/k$ times and 
${\sum }_i \ell _i = \mathcal {L}$ with the assumption of mutual independence of the 
$\ell _i$. Then, approximating discrete sums as integrals for large enough linking numbers, 
$\varOmega (\mathcal {L})$ is given by the volume of the simplex
\begin{equation} \varOmega(\mathcal{L}) = \int_{\ell_{min}}^{\mathcal{L}} \,\mathrm{d}\ell_1 \int_{\ell_{min}}^{\mathcal{L}} \,\mathrm{d}\ell_2 \cdots \int_{\ell_{min}}^{\mathcal{L}} \,\mathrm{d}\ell_k \,\delta\left( \sum_{i=1}^k \ell_i - \mathcal{L}\right) = \frac{(\mathcal{L}-\mathcal{L}_{min})^{k-1}}{(k-1)!}, \end{equation}
where the integration bounds correspond to the case of a positive helicity flow, in which linking numbers are bounded from below. Combining (4.4) and (4.5) while eliminating 
$\chi = (\langle \mathcal {L} \rangle - \mathcal {L}_{min})/k$ using (4.3), we obtain the k-Gamma distribution:
\begin{equation} p(\mathcal{L}) = \frac{k^k}{\varGamma(k)} \frac{\left( \mathcal{L} - \mathcal{L}_{min} \right)^{k-1}}{\left(\langle \mathcal{L} \rangle- \mathcal{L}_{min} \right)^k} \exp{\left(-k \frac{\mathcal{L} - \mathcal{L}_{min}}{\langle \mathcal{L} \rangle- \mathcal{L}_{min}}\right)}, \end{equation}
where 
$\varGamma$ is the Euler Gamma function. The k-Gamma distribution can be understood as a Gamma distribution with shape parameter 
$k$ for the scaled and shifted random variable 
$(\mathcal {L} - \mathcal {L}_{min})/(\langle \mathcal {L} \rangle - \mathcal {L}_{min})$. Note that in the case of negative helicity flows, we still obtain a k-Gamma distribution by substituting 
$\mathcal {L} \leftarrow - \mathcal {L}$ and 
$\mathcal {L}_{min} \leftarrow - \mathcal {L}_{max}$.
To test the validity of this approach, we fit the probability distribution function in (4.6) to the linking number distribution obtained from 
$N=1000$ vortex lines sampled from a positive helicity flow using maximum likelihood estimates of 
$\mathcal {L}_{min}$, 
$\langle \mathcal {L}\rangle$ and 
$k$ as a function of the vortex line length 
$L$ (figure 8). We sample 
$2\times 10^4$ points from the fitted distribution and find excellent agreement with the k-Gamma distribution, with 
$p(\mathcal {L})$ for long (
$L=100R$) vortex lines matching the fit down to shot noise. This increase in fit quality with increasing 
$L$ is consistent with the expectation that our independence and discrete sum approximations become more justified for longer vortex lines (figure 8a,b).
Figure 8. The linking number distribution at a fixed time is well described by a k-Gamma distribution. (a) For a flow with positive helicity (
$\varLambda = R/8$, 
$\kappa = 4/R$), the empirical probability distribution function of 
$\mathcal {L}$ (blue curves) for different vortex line lengths is well fit by a k-Gamma distribution (red curves), including up to shot noise for long enough vortex lines (log scale inset). Distributions are shown for 
$L/R=30, 40, 50, 60, 80, 100$, with 
$N=1000$ vortex lines. (b) Plotting centred and scaled distributions for 
$L/R \geq 50$ highlights the deviation from the normal distribution (black line) and shows the quality of the fit for large 
$L$. (c) The linking distribution has a Fano factor 
$\sigma ^2/\langle \mathcal {L} \rangle > 1$, showing super-Poissonian behaviour that is well captured by the k-Gamma fits in panel (c i). The skewness of fits asymptotically matches the data, which scales as 
$1/\sqrt {L}$, in agreement with the hypothesis of independent increments, shown in panel (c ii). (d) The fitted exponent 
$k$ scales linearly with 
$L$ (dashed lines as visual guide).
As expected from (3.5) and (4.1), the fitted value of 
$\langle \mathcal {L}\rangle$ recovers Arnold's equality. The fitted k-Gamma distributions also recover the correct scaling of mean and variance of data, notably displaying the super-Poissonian behaviour of the linking distribution as shown by its Fano factor 
$\sigma ^2/\langle \mathcal {L}\rangle >2$, and show the asymptotic behaviour of the skewness (figure 8c).
In previous applications of the k-Gamma distribution, 
$k$ is fit as a shape parameter determined by the variance 
$\sigma ^2$ of the distribution as 
$k = (\langle \mathcal {L}\rangle - \mathcal {L}_{min})/\sigma ^2$ (Aste et al. Reference Aste, Di Matteo, Saadatfar, Senden, Schröter and Swinney2007; Aste & Di Matteo Reference Aste and Di Matteo2008; Atia et al. Reference Atia2018; Day et al. Reference Day2022). To obtain parameter estimates as a function of vortex line length, we similarly fit our observed distribution using a maximum likelihood optimization. This procedure is agnostic to our arguments leading to (4.5), which posits that 
$k$ is set by the number of independent sub-domains of vortex lines. In support of the validity of the decomposition of vortex line linking numbers into contributions from sub-loops, we find that the fitted parameter 
$k$ linearly increases with the vortex line length 
$L$ (figure 8d). As we find that 
$k \approx L/(2R)$, we can estimate that topologically correlated domains have a characteristic size 
$\delta L \approx 2R$, indicating that this emergent correlation length is set by the diameter of the ball.
As we consider longer vortex lines, since 
$k\propto L$, we predict that the linking number distribution must eventually tend to the normal distribution. With 
$k\rightarrow +\infty$, the k-Gamma distribution converges to the Gaussian 
$\mathcal {N}(\mu =\langle \mathcal {L}\rangle, \sigma ^2)$, consistent with figure 8(b). This convergence is also expected of a sum of linking numbers 
$\ell _i$ from an increasingly large number of statistically independent sub-loops. For 
$\ell _i$ independent and identically distributed, the skewness of the resulting sum is expected to scale as 
$L^{-1/2}$ when each 
$\ell _i$ is drawn from a skewed distribution (Hall Reference Hall1992). Further validating the mesoscopic picture that we used to compute 
$\varOmega (\mathcal {L})$, this inverse-square root scaling behaviour is observed in figure 8(c).
5. Conclusion
Building on recent progress in spectral simulation techniques for spherical domains (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), we simulated a generalized Navier–Stokes (GNS) model for actively driven fluid flow in a confined 3-D domain to characterize the topology of the spontaneously forming chiral flow. Driven by a generic linear instability (Rothman Reference Rothman1989; Beresnev & Nikolaevskiy Reference Beresnev and Nikolaevskiy1993; Tribelsky & Tsuboi Reference Tribelsky and Tsuboi1996), we find that the GNS system produces an active Beltrami flow in the bulk for narrow energy injection bandwidth, a regime that could likely be realized in microfluidic experiments using semi-dense bacterial (Wioland et al. Reference Wioland, Woodhouse, Dunkel and Goldstein2016) or other microbial suspensions.
Leveraging recently published fast algorithms (Qu & James Reference Qu and James2021), we characterized the topological structure of this spontaneous flow through the pair-wise linking numbers of sampled vortex lines. We explicitly measured the convergence of the mean vortex line entanglement to the total helicity of the flow, as described asymptotically by Arnold. Importantly, these results apply to all simply confined solenoidal vector fields with appropriate boundary conditions, making our methodology applicable well beyond active matter models, including high-Reynolds-number and magnetohydrodynamic flows. In the active flow considered here, we found that a k-Gamma density describes the linking number distribution well, and propose a maximum entropy argument to explain this result.
We note that there has been substantial work on understanding the asymptotic behaviour of the winding number of two-dimensional random walks with or without chiral drift, in the presence of repulsion or confinement (Spitzer Reference Spitzer1958; Drossel & Kardar Reference Drossel and Kardar1996; Wen & Thiffeault Reference Wen and Thiffeault2019). In particular, it is known that the winding number of 2-D confined random walks eventually becomes normally distributed for long-enough walks. This convergence holds even in the presence of chiral drift, although the presence of drift leads to long-lasting skewness in the winding distribution (Drossel & Kardar Reference Drossel and Kardar1996; Wen & Thiffeault Reference Wen and Thiffeault2019). In comparison, the study of the linking number of confined random Brownian walks still presents many open questions (Orlandini & Whittington Reference Orlandini and Whittington2007), with to our knowledge no proof of convergence to normality, although results exist for the scaling of moments of drift-free polygonal Brownian walks (Panagiotou, Millett & Lambropoulou Reference Panagiotou, Millett and Lambropoulou2010; Marko Reference Marko2011). The numerically observed convergence to the normal distribution for the linking number of long confined vortex lines (figure 8) suggests connections between winding statistics and the linking numbers of random trajectories. Moreover, given the generic nature of our entropic reasoning, it is possible that the k-Gamma distribution is a universal feature of vortex line entanglement in chiral systems.
Finally, our analysis shows that even a small number of vortex lines with lengths 
$L \sim R$ are sufficient to infer both the helicity sign and the helicity value. This rapid convergence indicates that observations of vortex line entanglement – for instance through tracer particles (Kleckner & Irvine Reference Kleckner and Irvine2013) or embedded filaments (Kirchenbuechler et al. Reference Kirchenbuechler, Guu, Kurniawan, Koenderink and Lettinga2014; Du Roure et al. Reference Du Roure, Lindner, Nazockdast and Shelley2019) – could provide reliable measurements of chiral symmetry breaking in 3-D experimental flows even with limited data.
Acknowledgements
The authors thank M. Kardar for valuable discussions, and the MIT SuperCloud and Lincoln Laboratory Supercomputing Center for providing HPC resources that have contributed to the results reported within this paper.
Funding
This work was supported by a MathWorks Science Fellowship (N.R.), NSF Award DMS-1952706 (N.R. and J.D.), Alfred P. Sloan Foundation Award G-2021-16758 (J.D.), the Robert E. Collins Distinguished Scholarship Fund (J.D.) and the MathWorks Professorship Fund (J.D.) of the MIT Mathematics Department, the Swiss National Science Foundation Ambizione Grant PZ00P2_202188 (J.S.), and NASA HTMS Award 80NSSC20K1280 (K.J.B.).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The simulation and linking analysis scripts are available at the following repositories: Simulation – https://github.com/kburns/gns_ball; Analysis – https://github.com/NicoRomeo/pyfln.
Appendix. Numerical methods
We solve the GNS equations (1.3) and (2.2) using Dedalus v3 (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). The variables are discretized using the full-ball basis from Vasil et al. (Reference Vasil, Lecoanet, Burns, Oishi and Brown2019); Lecoanet et al. (Reference Lecoanet, Vasil, Burns, Brown and Oishi2019), consisting of spin-weighted spherical harmonics in the angular directions and one-sided Jacobi polynomials in radius. This discretization satisfies the exact regularity conditions of smooth tensor fields at the poles and the origin, without having to excise any regions around the coordinate singularities. Other attempts to simulate the GNS equations in the ball using alternate discretizations have exhibited stability issues (Boullé, Słomka & Townsend Reference Boullé, Słomka and Townsend2021), which may be related to their inexact imposition of regularity conditions at the origin.
In Dedalus, the equations are integrated semi-implicitly with the linear terms solved implicitly and the nonlinear computed explicitly using a third-order four-stage mixed Runge–Kutta (RK) scheme (Ascher et al. Reference Ascher, Ruuth and Spiteri1997, § 2.8). The nonlinear terms are computed pseudospectrally with 3/2 dealiasing. Incompressibility is directly enforced by implicitly solving for the dynamic pressure as a Lagrange multiplier in each RK substep. The boundary conditions are implemented using a generalized tau method, where radial tau terms are explicitly added to the linear system. The system is evolved with a fixed time step from random initial noise until reaching a statistically stationary state.
The initial conditions for velocity are Gaussian white fields with mean value 
$\langle \boldsymbol {u}\rangle = 0$ and spatial correlations between components of the velocity 
$\langle u_i(\boldsymbol {x}) u_j(\kern0.7pt \boldsymbol {y})\rangle = \sigma ^2 \delta _{ij}\delta (\boldsymbol {x}-\boldsymbol {y})$, where 
$i,j$ index the velocity components. The random field is then projected onto its incompressible component such that 
$\boldsymbol {\nabla } \boldsymbol{\cdot} \boldsymbol {u} = 0$. We find that the statistical steady state is independent of the amplitude of the initial amplitude 
$\sigma$ (figure 9), consistent with the fact that the initial transient regime is dominated by a linear instability.
Figure 9. The flow statistics in steady state are independent of initial condition amplitude. Energy and total helicity as a function of time for random initial conditions with different amplitudes 
$\sigma$. All simulations have 
$\varLambda = R/8, \kappa =4/R, \tau = 1$. (a) 20 simulations with 
$\sigma = 10^{-2}$. (b) 20 simulations with 
$\sigma = 10^{-4}$. The bottom panel shows the energy on a semi-logarithmic scale, indicating an exponential growth of the energy at early times, characteristic of a linear instability. (c) 20 simulations with 
$\sigma = 10^{-5}$. We note that the main text results are all reported for 
$\sigma = 10^{-3}$. All values are in simulation units where 
$R= 1, \tau = 1$.
Vortex line integration is performed using Scipy's ordinary differential equation integrator called through the method 
$\texttt {solve_ivp}$. The RK45 solver is used with absolute tolerance 
$10^{-6}$, relative tolerance 
$10^{-12}$ and maximal step size 
$10^{-2}$ (in length units where the ball radius is 
$R=1$). Linear interpolation of the vorticity field between points on the numerical quadrature grid is used to integrate the vortex lines.
