2.1 Conservation equations for active nematohydrodynamics

In this formulation, the mesoscopic orientational order of active particles is represented by the nematic tensor $\unicode[STIX]{x1D618}_{ij}=(3q/2)(n_{i}n_{j}-\unicode[STIX]{x1D6FF}_{ij}/3)$ , where $q$ is the magnitude of the nematic order, $n_{i}$ is the director and $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta. The conservation equation for the nematic tensor is given by

(2.1) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D618}_{ij}+u_{k}\unicode[STIX]{x2202}_{x_{k}}\unicode[STIX]{x1D618}_{ij}=\unicode[STIX]{x1D6E4}\unicode[STIX]{x1D60F}_{ij}+\unicode[STIX]{x1D619}_{ij},\end{eqnarray}$$

where $u_{k}$ are velocity components, $t$ is time, $x_{i}$ are spatial coordinates and $\unicode[STIX]{x1D6E4}$ is proportional to a rotational diffusivity. Additionally, $\unicode[STIX]{x1D619}_{ij}$ is a standard co-rotation tensor defined as

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D619}_{ij}=(\unicode[STIX]{x1D706}\unicode[STIX]{x1D61A}_{ik}+\unicode[STIX]{x1D6FA}_{ik})\left(\unicode[STIX]{x1D618}_{kj}+\frac{\unicode[STIX]{x1D6FF}_{kj}}{3}\right)+\left(\unicode[STIX]{x1D618}_{ik}+\frac{\unicode[STIX]{x1D6FF}_{ik}}{3}\right)(\unicode[STIX]{x1D706}S_{kj}-\unicode[STIX]{x1D6FA}_{kj})-2\unicode[STIX]{x1D706}\left(\unicode[STIX]{x1D618}_{ij}+\frac{\unicode[STIX]{x1D6FF}_{ij}}{3}\right)\unicode[STIX]{x1D618}_{kl}\unicode[STIX]{x2202}_{x_{k}}u_{l},\end{eqnarray}$$

which accounts for the response of the orientation field to the extensional and rotational components of the velocity gradients, with $\unicode[STIX]{x1D61A}_{ij}=(1/2)(\unicode[STIX]{x2202}_{x_{i}}u_{j}+\unicode[STIX]{x2202}_{x_{j}}u_{i})$ and $\unicode[STIX]{x1D6FA}_{ij}=(1/2)(\unicode[STIX]{x2202}_{x_{j}}u_{i}-\unicode[STIX]{x2202}_{x_{i}}u_{j})$ being the strain-rate and vorticity tensors, respectively. The relative importance of the vorticity and strain rate is controlled by the parameter $\unicode[STIX]{x1D706}$ , which characterizes the alignment of the nematics with the flow. The term involving $\unicode[STIX]{x1D6E4}$ in (2.1) describes, through the auxiliary tensor  the relaxation of $\unicode[STIX]{x1D618}_{ij}$ to a minimum of a free energy

(2.3) $$\begin{eqnarray}{\mathcal{F}}=(A/2)\unicode[STIX]{x1D618}_{ij}\unicode[STIX]{x1D618}_{ij}+(B/3)\unicode[STIX]{x1D618}_{ij}\unicode[STIX]{x1D618}_{jk}\unicode[STIX]{x1D618}_{ki}+(C/4)(\unicode[STIX]{x1D618}_{ij}\unicode[STIX]{x1D618}_{ij})^{2}+(K/2)(\unicode[STIX]{x2202}_{x_{k}}\unicode[STIX]{x1D618}_{ij})^{2}.\end{eqnarray}$$

In this formulation, represents the variational derivative, $\text{tr}$ denotes the trace, $K$ is an elastic constant, while $A$ , $B$ and $C$ are material constants that determine the equilibrium state of the orientational order. The first three terms in (2.3) correspond to the Landau/De Gennes free energy, while the last term represents the Frank elastic energy with a one-constant approximation (De Gennes & Prost Reference De Gennes and Prost1995). Equation (2.3) is the free energy expansion in terms of the order parameter, where the terms allowed by symmetry have been retained. Note that the Frank term in (2.3) gives rise to a diffusive flux of $\unicode[STIX]{x1D618}_{ij}$ in the form $\unicode[STIX]{x1D6E4}K\unicode[STIX]{x2202}_{x_{k},x_{k}}^{2}\unicode[STIX]{x1D618}_{ij}$ in (2.1). The description of the flow field is completed by the mass and momentum conservation equations, namely

(2.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{x_{i}}u_{i}=0,\quad \unicode[STIX]{x1D70C}\unicode[STIX]{x2202}_{t}u_{i}+\unicode[STIX]{x1D70C}u_{j}\unicode[STIX]{x2202}_{x_{j}}u_{i}=-\unicode[STIX]{x2202}_{x_{i}}p+\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{x_{j},x_{j}}^{2}u_{i}+\unicode[STIX]{x2202}_{x_{j}}\unicode[STIX]{x1D70E}_{ij}-\unicode[STIX]{x1D701}\unicode[STIX]{x2202}_{x_{j}}\unicode[STIX]{x1D618}_{ij},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}$ is the density, $p$ is the pressure and $\unicode[STIX]{x1D707}$ is the dynamic viscosity. The additional terms in the momentum equation (2.4) involve the elastic stress

(2.5)

which represents the passive conformational resistance of the nematics to changes in the orientational order (Edwards, Beris & Grmela Reference Edwards, Beris and Grmela1991), and the active stress $\unicode[STIX]{x1D701}\unicode[STIX]{x1D618}_{ij}$ , which corresponds to a coarse-grained collective effect of the stresslets set up by the active particles (Simha & Ramaswamy Reference Simha and Ramaswamy2002). In particular, the divergence of the active stress is responsible for the injection of kinetic energy through the orientational order of particles represented by the nematic tensor, whose conservation equation (2.1) is nonlinearly coupled with the velocity field. The absolute value of the coefficient $\unicode[STIX]{x1D701}$ is proportional to the activity, where positive and negative values of $\unicode[STIX]{x1D701}$ correspond, respectively, to extensile (pusher) and contractile (puller) particles, the former being the case addressed in the computations shown below.

2.2 Remarks on the conservation equations

Some aspects with regard to the formulation given above are worth briefly pointing out. The first one is related to characteristic dimensionless parameters. In the conditions addressed by these simulations (see further details in § 2.3) and most experimental observations, the Reynolds number is small, $Re=\unicode[STIX]{x1D70C}u^{\prime }\ell /\unicode[STIX]{x1D707}\lesssim 1$ , with $u^{\prime }$ and $\ell$ being the integral scales for velocity and length, respectively. In this limit, the convective transport of momentum has a diminishing influence, with the dominant balance in (2.4) being between viscous and active stresses. Similarly, the orientational order diffuses slower than momentum, as indicated by the relatively high Schmidt number $Sc=(\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C})/(\unicode[STIX]{x1D6E4}K)\gtrsim 5$ , which results in characteristic structures of the velocity field that are larger than those of the orientational order, as discussed in § 3. As a consequence, the Péclet number of the nematic order $Pe=ReSc\gtrsim 3$ suggests that advection may play a more important role in transporting orientational order than it does for transporting momentum.

The second aspect worth stressing is that the conservation equations outlined above do not include any obvious external forcing term aimed at sustaining the dynamics. A flow laden with active matter is different from an inactive flow externally driven by boundary conditions, imposed shear or volumetric forces added to the momentum conservation equation. In real biological flows laden with swimming bacteria, adenosine tri-phosphate (ATP) molecules are consumed by the bacteria via hydrolysis reactions and transformed into motion in such a way that the driving occurs internally at the expense of depletion of chemical energy. The source of chemical energy, however, is not represented in the formulation above, in that it only concerns systems where the rate of ATP depletion is infinitesimally small compared to hydrodynamic processes.

The active motorization of the flow can be understood by multiplying the momentum equation in (2.4) by $u_{i}$ and performing a surface (in two dimensions) or volumetric (in three dimensions) periodic integration, which leads to the balance equation $\unicode[STIX]{x1D70C}(\text{d}k/\text{d}t)=-\unicode[STIX]{x1D716}-\langle \unicode[STIX]{x1D70E}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle +\unicode[STIX]{x1D701}\langle \unicode[STIX]{x1D618}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle$ for the spatially averaged kinetic energy $k=\langle u_{i}u_{i}/2\rangle$ , where $\unicode[STIX]{x1D716}=\langle 2\unicode[STIX]{x1D707}\unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle$ is the viscous dissipation. The work done by the active forces, given by the term $\unicode[STIX]{x1D701}\langle \unicode[STIX]{x1D618}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle$ , represents the main source of $k$ . The power deployed by the active work is dissipated by viscosity as shown below, thereby yielding a stationary state in which $k$ remains mostly constant.

Figure 1. Instantaneous contours of (a) velocity magnitude (high activity, main panel; low activity, inset), (b) vorticity and (c) magnitude of the nematic tensor. In (b,c), which are zoomed views of the white squared region in (a), green lines are nematic director fields while symbols represent $+1/2$ (circles) and $-1/2$ (triangles) topological defects. The small inset above (c) shows the director field around topological defects.

2.3 Computational set-up

The formulation described above is integrated numerically in quasi-two-dimensional (quasi-2-D) and 3-D domains with periodic boundary conditions. In the quasi-2-D framework the velocity varies in two dimensions while the nematic directors could develop out-of-plane components, as in previous studies (Thampi, Golestanian & Yeomans Reference Thampi, Golestanian and Yeomans2013; Doostmohammadi et al. Reference Doostmohammadi, Adamer, Thampi and Yeomans2016b ; Saw et al. Reference Saw, Doostmohammadi, Nier, Kocgozlu, Thampi, Toyama, Marcq, Lim, Yeomans and Ladoux2017). The dimensional parameters in the simulations are $A=0.04$ , $B=0$ for two dimensions and $B=0.06$ for three dimensions, $C=0.06$ , $\unicode[STIX]{x1D706}=0.3$ (tumbling particles), $\unicode[STIX]{x1D6E4}=0.34$ , $K=0.40$ , $\unicode[STIX]{x1D70C}=1.0$ and $\unicode[STIX]{x1D707}=0.66$ ( $Sc=4.90$ ). The baseline activity coefficient is $\unicode[STIX]{x1D701}_{0}=0.036$ for two and three dimensions, with an additional 2-D computation being performed with a smaller activity $\unicode[STIX]{x1D701}_{0}/10$ . All parameters here are expressed in simulation units. This generic set of parameters has been shown in earlier work to reproduce qualitatively the active turbulent state observed in microtubule/motor-protein suspensions (Sanchez et al. Reference Sanchez, Chen, DeCamp, Heymann and Dogic2012; Wensink et al. Reference Wensink, Dunkel, Heidenreich, Drescher, Goldstein, Lowen and Yeomans2012; Thampi et al. Reference Thampi, Golestanian and Yeomans2013). A standard lattice-Boltzmann approach is used to integrate (2.4), while (2.1) is solved by employing a second-order finite-difference predictor–corrector algorithm. The resulting set of ordinary differential equations is integrated in time using an Euler method. The number of grid points for the 2-D and 3-D simulations is $N=512^{2}$ and $128^{3}$ , respectively, with a minimum grid spacing of $\unicode[STIX]{x1D6E5}=\ell _{d}/3$ for all cases, where $\ell _{d}=\sqrt{K/A}=3.16$ is the characteristic size of the topological defects in the orientational order field. The time increment used in the numerical integrations is $\unicode[STIX]{x0394}t=t_{q}/60$ , with $t_{q}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D701}$ being a characteristic time scale for the damping of the activity by viscosity. The initial conditions consist of zero velocity everywhere while the directors are set to random orientations. Data are collected after approximately $2\times 10^{4}$ time steps for 10 snapshots spanning a time period of $150t_{q}$ , which provide $512^{2}\times 10$ and $128^{3}\times 10$ statistical samples leading to converged probability density functions (PDFs) in two and three dimensions, respectively. In the results, spatial coordinates are normalized with $\ell _{d}$ . Additionally, velocities are normalized with $u^{\prime }=0.11$ (for two dimensions with $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{0}$ ) or $u^{\prime }=0.04$ (for three dimensions, and for two dimensions with $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{0}/10$ ), where $u^{\prime }=\sqrt{2k}$ .

3.1 Flow structures

Instantaneous contours of velocity and vorticity are provided in figures 1(a,b) and 2(a) for 2-D and 3-D domains, respectively. Specifically, the effects of decreasing the activity coefficient from $\unicode[STIX]{x1D701}_{0}$ to $\unicode[STIX]{x1D701}_{0}/10$ in the 2-D simulations are an increase in the integral length (computed from the ensemble-averaged kinetic-energy spectrum) from $\ell =5.45$ to $\ell =10.22$ , a decrease in the Reynolds number from to $Re=0.74$ to $Re=0.68$ and a decrease in the dissipation from $\unicode[STIX]{x1D716}=1\times 10^{-3}$ to $\unicode[STIX]{x1D716}=2\times 10^{-4}$ . The resulting low-activity flow has a less dense pattern of flow structures, as observed by comparing the main (high activity) and inset (low activity) frames in figure 1(a) and noticing that the normalizing length $\ell _{d}$ is independent of activity. It is noteworthy that, for the same activity coefficient $\unicode[STIX]{x1D701}_{0}$ , moving from two to three dimensions translates into an increase in the integral length from $\ell =5.45$ to $\ell =8.97$ , a decrease in the Reynolds number from $Re=0.74$ to $Re=0.53$ and a decrease in the dissipation from $\unicode[STIX]{x1D716}=1\times 10^{-3}$ to $\unicode[STIX]{x1D716}=0.3\times 10^{-4}$ .

Figure 2. (a) Instantaneous isosurfaces of velocity magnitude $(u_{i}u_{i})^{1/2}/u^{\prime }=1.1$ (30 % of maximum value) and enstrophy $\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{i}(\ell _{d}/u^{\prime })^{2}=0.4$ (20 % of maximum value) for 3-D simulations. (b) Scatter plot of the velocity-gradient invariants. (c) Schematics of the stresslet-like flow induced by extensile active particles ( $\unicode[STIX]{x1D701}>0$ ).

The spatial variations in the nematic tensor $\unicode[STIX]{x1D618}_{ij}$ are central to the generation of vorticity. This is easily observed by taking the curl of (2.4), namely

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D70C}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D714}_{i}+\unicode[STIX]{x1D70C}u_{j}\unicode[STIX]{x2202}_{x_{j}}\unicode[STIX]{x1D714}_{i}=-\unicode[STIX]{x1D70C}\unicode[STIX]{x1D714}_{j}\unicode[STIX]{x2202}_{x_{j}}u_{i}+\unicode[STIX]{x1D700}_{ijk}\unicode[STIX]{x2202}_{x_{j}}(\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{x_{\ell },x_{\ell }}^{2}u_{k}+\unicode[STIX]{x2202}_{x_{\ell }}\unicode[STIX]{x1D70E}_{k\ell }-\unicode[STIX]{x1D701}\unicode[STIX]{x2202}_{x_{\ell }}Q_{k\ell }),\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{i}$ is the vorticity and $\unicode[STIX]{x1D700}_{ijk}$ is the permutation tensor. In two dimensions, the vortex stretch term is exactly zero and the dominant mechanism of vorticity generation is the curl of the divergence of the active stresses. The structures of vorticity, which are shown in figure 1(b), are different from the classic round vortices observed in high-Reynolds-number 2-D isotropic flows. Instead, vortical structures attain here band-like shapes, which are closely related to thinner elongated regions referred to as walls, where the magnitude of the nematic-order tensor becomes small, as shown by the solid contours of $\unicode[STIX]{x1D714}_{3}$ overlaid on the director field (largest eigenvector of the $\unicode[STIX]{x1D618}_{ij}$ tensor) in figure 1(b). The walls are characterized by bend deformations in the in-plane director field in figure 1(c) (note that the out-of-plane components in these 2-D simulations are zero), which separate nematically aligned regions ( $q\sim 1$ ) across interstitial isotropic states ( $q\ll 1$ ), and are typically much thinner than the hydrodynamic structures of velocity and vorticity. The resulting $\pm 1/2$ topological defects, depicted by circles and triangles in figure 1(b,c) (see inset), represent singular, disordered regions of strong vorticity generation that are created and annihilate in pairs while propagating in the flow in a complex manner (Doostmohammadi et al. Reference Doostmohammadi, Thampi and Yeomans2016a ; Saw et al. Reference Saw, Doostmohammadi, Nier, Kocgozlu, Thampi, Toyama, Marcq, Lim, Yeomans and Ladoux2017) that is beyond the scope of the present study.

In three dimensions, the vortex stretch term in (3.1) represents a much smaller contribution than the active stresses because of the low Reynolds numbers involved. As observed in figure 2(a), the resulting vortical structures in three dimensions are elongated as well but smaller than the velocity ones. The 3-D mechanisms of vorticity generation are mostly unknown since the description of 3-D topological defects in active nematics is not well understood. Further insight into rotational and straining components of the 3-D flow field can be gained by examining the velocity-gradient invariants $Q_{inv}=(1/4)(\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{i}-2\unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij})$ and $R_{inv}=(3/4)(\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{j}\unicode[STIX]{x1D61A}_{ij}+4\unicode[STIX]{x1D61A}_{ij}S_{jk}S_{ki})$ , which are expedient for classifying flow structures. In contrast to typical scatter plots for high-Reynolds-number turbulence, where most of the activity is in the upper-right and lower-left quadrants, figure 2(b) shows that for active flows straining is predominant. As a result, the marginal PDF of $Q_{inv}$ is heavily skewed toward negative values (skewness $-1.6$ ). The straining is caused by the cumulative effect of the stresslets from the extensile active particles (see flow sketch in figure 2 c).

Figure 3. Ensemble-averaged PDFs of (a) velocity, (b) velocity gradient, (c) vorticity and (d) nematic-order magnitude, including three dimensions ( $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{0}$ , solid lines) and two dimensions ( $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{0}$ , dot-dashed lines; $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{0}/10$ , dashed lines). Red short dashed lines indicate reference Gaussian distributions.

3.2 PDF moments of flow variables

The PDFs of velocity $u_{1}$ , velocity gradient $\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}$ , vorticity $\unicode[STIX]{x1D714}_{3}$ and magnitude of the nematic order $q$ are shown in figure 3, and some of their moments are listed in table 1. The PDFs of velocity and vorticity have nearly Gaussian flatness, while the largest skewness of the velocity gradient is reached in the 2-D high-activity case and equals $-0.10$ . In all cases, the vorticity flatness and the velocity-gradient skewness are smaller than typical values observed in high-Reynolds-number turbulence (i.e.  ${\sim}8$ and ${\sim}\!-0.4$ , respectively). In two dimensions, small activities favour sub-Gaussian flatness for vorticity and velocity along with decreasing skewness of the velocity gradient. Note, however, that for the same activity coefficient the 3-D case leads to a smaller skewness of the velocity gradient, which suggests that the lesser spatial confinement plays a role in the development of fluctuations. Additionally, small values of the nematic-order magnitude, which correspond to isotropic behaviour and vorticity generation, are statistically favoured at large activities.

Table 1. Moments of velocity, velocity derivative and vorticity PDFs.

3.3 Intermittency analysis

Albeit small, intermittency may not be entirely ruled out in these systems, as suggested by a narrowband filtering analysis of the results based on a discrete db-4 wavelet decomposition of the velocity and vorticity fields. This is shown in figure 4, which provides the scale-dependent flatness $F^{(s)}$ of the PDFs of the direction-averaged velocity and vorticity wavelet coefficients, $\check{u}_{1}^{(s)}$ and $\check{\unicode[STIX]{x1D714}}_{3}^{(s)}$ , normalized by their corresponding standard deviations $\unicode[STIX]{x1D70E}_{u}^{(s)}$ and $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}}^{(s)}$ , with $s$ being a scale index that ranges from $1$ (equivalent to a length scale of $2\unicode[STIX]{x1D6E5}$ ) to $s_{max}=\log _{2}N$ (equivalent to $N\unicode[STIX]{x1D6E5}$ ) and is related to a representative wavenumber as $\unicode[STIX]{x1D705}=2\unicode[STIX]{x03C0}2^{-s}/\unicode[STIX]{x1D6E5}$ . In these simulations, $s_{max}=9$ (for 2-D cases) and $s_{max}=7$ (for 3-D cases). The inverted caret denotes the wavelet transform, e.g.  $\check{u}_{i}^{(s,d)}(\boldsymbol{x}_{\boldsymbol{s}})=\langle u_{i}(\boldsymbol{x})\unicode[STIX]{x1D6F9}^{(s,d)}(\boldsymbol{x}-\boldsymbol{x}_{\boldsymbol{s}})\rangle$ , where $d$ is a positive-integer direction index ( $d=1,2,3$ and $d=1,2,\ldots ,7$ in two and three dimensions, respectively), $\boldsymbol{x}_{\boldsymbol{s}}=2^{s-1}(i\unicode[STIX]{x1D6E5},j\unicode[STIX]{x1D6E5},k\unicode[STIX]{x1D6E5})$ are scale-dependent wavelet grids where the wavelets are collocated, with $i$ , $j$ , $k=1,3,5,\ldots ,N/2^{s-1}-1$ , and $\unicode[STIX]{x1D6F9}^{(s,d)}(\boldsymbol{x}-\boldsymbol{x}_{\boldsymbol{s}})$ are wavelet basis functions that are here taken to be tensor products of orthonormal 1-D db-4 wavelets with four vanishing moments. The reader is referred to Meneveau (Reference Meneveau1991) and Schneider & Vasilyev (Reference Schneider and Vasilyev2010) for general applications of wavelets in turbulent flows, and to Nguyen et al. (Reference Nguyen, Farge and Schneider2012) for a scale-dependent flatness analysis of homogeneous isotropic turbulence similar to the one performed here.

While the 2-D fields remain nearly Gaussian at all scales, with a slight increase for the vorticity flatness observed in the largest scale, figure 4(b) shows that the 3-D fields contain significant intermittency in the small scales, as indicated by the strong increase in $F_{s}$ with $\unicode[STIX]{x1D705}$ (main panel) and by the increasingly longer tails in the PDFs of the wavelet coefficients as the length scale increases (inset). Nonetheless, the energetic content of these small scales and their associated intermittent motion is small in all cases. This can be understood by noticing the rapid decay of the Fourier kinetic-energy and enstrophy spectra, $E_{k}$ and $E_{\unicode[STIX]{x1D714}}$ , as $\unicode[STIX]{x1D705}$ increases, as observed in figure 5(a,b,d,e). Specifically, the kinetic-energy spectra decay with a slope that decreases from $-4.5$ to $-3.5$ as the activity is increased tenfold. As a result, the enstrophy spectra reach a maximum at the integral scale and decay rapidly thereafter, indicating that the small-scale gradients bear vanishing energy. This detracts dynamical relevance from the increased intermittency observed in the 3-D small scales and sets a fundamental difference between these flows and high-Reynolds-number turbulence; in the latter, the kinetic-energy spectra decays at a slower rate and the small-scale vorticity intermittency is energetic.

It is also of interest to note that the characteristic wavenumber where the nematic-order fluctuation energy spectra $E_{q}$ (defined such that the area under the curve is $\langle q^{\prime }q^{\prime }\rangle$ ) reach a maximum is approximately one decade smaller than the wavenumber corresponding to the maximum enstrophy spectra in the 2-D low-activity case, as shown in figure 5(c). The distance between peaks decreases with increasing activity and moving to three dimensions, as observed in figure 5(f). The wavenumber of maximum $E_{q}$ decreases as the activity increases and its value is closer to $2\unicode[STIX]{x03C0}/\ell _{d}$ than to $2\unicode[STIX]{x03C0}/\ell$ , which spectrally illustrates the fine structure of the nematic-order field compared to the coarser velocity field.

Figure 4. Wavelet-based scale-dependent flatness (main frames) for velocity and vorticity in (a) 2-D and (b) 3-D cases at $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{0}$ , along with corresponding scale-conditioned PDFs (insets) of the wavelet coefficients. The PDFs are conditioned on $s=1$ (dark blue lines) $s=2$ (red), $s=3$ (orange), $s=4$ (purple), $s=5$ (green) and $s=6$ (light blue).

Figure 5. Ensemble-averaged Fourier and wavelet spectra of kinetic energy. The solid contours in (d–f) correspond to the PDF of the wavelet spectra, which include the mean (solid lines) and associated 95 % confidence intervals (dashed lines).

Figure 6. Mean and 95 % confidence intervals of wavelet-based spectral energy-transfer flux of the 3-D flow for (a) active and convective fluxes and (b) viscous flux. The panels also show solid contours for PDFs of (a) active and (b) viscous fluxes, indicating spatial variabilities.

3.4 Spectral energy-transfer analysis

As discussed in § 2.2, a crucial role in the dynamics is played by the divergence of the active stress $-\unicode[STIX]{x1D701}\unicode[STIX]{x1D618}_{ij}$ . Because of the small $Re$ involved, it is anticipated that the spectral transfer of the active energy is locally dissipated by viscosity, since the convective inter-scale transfer is a mechanism of secondary importance. Although there may exist additional triadic interactions resulting from (2.5) that could transport energy across scales, the 3-D results provided in figure 6 for active ( $\check{{\mathcal{T}}}_{A}$ ), convective ( $\check{{\mathcal{T}}}_{C}$ ) and viscous ( $\check{{\mathcal{T}}}_{V}$ ) wavelet-based spectral energy-transfer fluxes support the view that locality may dominate the transfer. Specifically, these fluxes describe the rate at which the spatially averaged, wavelet spectral kinetic-energy density, $E_{k}=2^{-3s}\langle \sum _{d}\check{u}_{i}^{(s,d)}(\boldsymbol{x}_{\boldsymbol{s}})\check{u}_{i}^{(s,d)}(\boldsymbol{x}_{\boldsymbol{s}})/2\rangle _{\boldsymbol{x}_{\boldsymbol{s}}}/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}$ , is transferred across scales. In particular, $E_{k}$ and the analogous $E_{\unicode[STIX]{x1D714}}$ and $E_{q}$ are shown and compared to their Fourier-based counterparts in figure 5(d–f). In this formulation, $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}=2\unicode[STIX]{x03C0}\ln 2/(2^{s}\unicode[STIX]{x1D6E5})$ is a discrete wavenumber shell, and the subindexed bracketed operator represents spatial averaging over scale-dependent wavelet grids $\boldsymbol{x}_{\boldsymbol{s}}$ .

Upon wavelet-transforming the momentum equation in (2.4), multiplying by $\check{u}_{i}^{(s,d)}(\boldsymbol{x}_{\boldsymbol{s}})$ and summing over $d$ , the spectral-energy equation $\unicode[STIX]{x2202}E_{k}/\unicode[STIX]{x2202}t=\sum \check{{\mathcal{T}}}(\unicode[STIX]{x1D705})$ is obtained. Here, the source term represents the sum of spectral energy-transfer fluxes created by each term on the right-hand side of the momentum equation in (2.4). For any force $\unicode[STIX]{x1D719}_{i}$ , the corresponding spectral flux is given by $\check{{\mathcal{T}}}(\unicode[STIX]{x1D705})=[(2^{-2s}\unicode[STIX]{x1D6E5})/(2\unicode[STIX]{x03C0}\ln 2)]\langle \sum _{d}\check{u}_{i}^{(s,d)}(\boldsymbol{x}_{\boldsymbol{s}})\check{\unicode[STIX]{x1D719}}_{i}^{(s,d)}(\boldsymbol{x}_{\boldsymbol{s}})\rangle _{\boldsymbol{x}_{\boldsymbol{s}}}$ , with $\check{{\mathcal{T}}}>0$ and $\check{{\mathcal{T}}}<0$ indicating, respectively, inflow and outflow of energy at a given $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}$ . Note that $\check{{\mathcal{T}}}=\check{{\mathcal{T}}}_{A}$ for $\unicode[STIX]{x1D719}_{i}=-\unicode[STIX]{x1D701}\unicode[STIX]{x2202}_{x_{j}}\unicode[STIX]{x1D618}_{ij}$ , $\check{{\mathcal{T}}}=\check{{\mathcal{T}}}_{V}$ for $\unicode[STIX]{x1D719}_{i}=\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{x_{j},x_{j}}^{2}u_{i}$ and $\check{{\mathcal{T}}}=\check{{\mathcal{T}}}_{C}$ for $\unicode[STIX]{x1D719}_{i}$ $=-\unicode[STIX]{x1D70C}u_{j}\unicode[STIX]{x2202}_{x_{j}}u_{i}$ , which satisfy $\sum _{\unicode[STIX]{x1D705}}\check{{\mathcal{T}}}_{A}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}=\unicode[STIX]{x1D701}\langle \unicode[STIX]{x1D618}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle$ , $\sum _{\unicode[STIX]{x1D705}}\check{{\mathcal{T}}}_{V}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}=-\unicode[STIX]{x1D716}$ and $\sum _{\unicode[STIX]{x1D705}}\check{{\mathcal{T}}}_{C}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}=0$ .

Figure 6(a) indicates that the active stress acts as a kinetic-energy source at all scales on spatial average, with the maximum mean of $\check{{\mathcal{T}}}_{A}$ occurring at scales of the same order as the integral length. Conversely, the viscous flux $\check{{\mathcal{T}}}_{V}$ is a sink of kinetic energy and has a trend that is exactly opposite to $\check{{\mathcal{T}}}_{A}$ , as shown in figure 6(b), which indicates that the active energy is mostly dissipated locally in spectral space by viscosity. The spatial localization of the transfer, which is illustrated by the unbracketed versions of $E_{k}$ , $E_{\unicode[STIX]{x1D714}}$ , $E_{q}$ and $\check{{\mathcal{T}}}$ , is represented by the variabilities of the PDFs shown in figures 5(d–f) and 6. Specifically, the PDFs in figure 6 reveal that the viscous transfer flux is spatially correlated with the active one (correlation coefficient $-0.73$ at $s=4$ ), suggesting that upon deployment the active energy is dissipated mostly locally also in physical space.

The physical picture implied by figure 6 provides no evidence for an energy cascade in momentum where the sink $\unicode[STIX]{x1D716}$ and main source $\unicode[STIX]{x1D701}\langle \unicode[STIX]{x1D618}_{ij}\unicode[STIX]{x1D61A}_{ij}\rangle$ of mean kinetic energy could act in disparate ranges of scales interacting through a crossing long-range mechanism. This is in contrast to high-Reynolds-number turbulence and its clear separation of scales between the large-scale forcing range and the small-scale molecular dissipation range. These conclusions could, however, be different for the nematic-order energy, in that the transport description of the latter is highly nonlinear and involves cross-triadic terms with the velocity as in (2.2). These aspects will be the subject of future research.