# An analytical and computational study of the incompressible Toner-Tu
Equations

John D. Gibbon<EMAIL_ADDRESS>Department of Mathematics, Imperial
College London, London SW7 2AZ, UK  Kolluru Venkata Kiran
<EMAIL_ADDRESS>Centre for Condensed Matter Theory, Department of
Physics, Indian Institute of Science, Bangalore, 560012, India  Nadia Bihari
Padhan<EMAIL_ADDRESS>Rahul Pandit<EMAIL_ADDRESS>
###### Abstract

The incompressible Toner-Tu (ITT) partial differential equations (PDEs) are an
important example of a set of active-fluid PDEs. While they share certain
properties with the Navier-Stokes equations (NSEs), such as the same scaling
invariance, there are also important differences. The NSEs are usually
considered in either the decaying or the additively forced cases, whereas the
ITT equations have no additive forcing. Instead, they include a linear,
activity term $\alpha\mbox{\boldmath$u$}$ ($u$ is the velocity field) which
pumps energy into the system, but also a negative
$\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2}$-term which provides a platform
for either frozen or statistically steady states. Taken together, these
differences make the ITT equations an intriguing candidate for study using a
combination of PDE analysis and pseudo-spectral direct numerical simulations
(DNSs). In the $d=2$ case, we have established global regularity of solutions,
but we have also shown the existence of bounded hierarchies of weighted, time-
averaged norms of both higher derivatives and higher moments of the velocity
field. Similar bounded hierarchies for Leray-type weak solutions have also
been established in the $d=3$ case. We present results for these norms from
our DNSs in both $d=2$ and $d=3$, and contrast them with their Navier-Stokes
counterparts.

###### keywords:

Active matter, Toner-Tu, Navier-Stokes, PDE

††journal: Physica D

##  1 Introduction

The field of active matter continues to grow rapidly [2, 3, 4, 5, 6, 7, 8, 9,
10, 11, 12, 14, 15, 16, 17, 18, 13, 19, 20, 21, 23, 22, 24, 25, 26]. The term
is generally used for systems that have bodies, e.g. birdoids in computer
animations [2], birds in a flock [3, 5, 6], or bacteria in dense suspensions
[5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 13, 16, 17, 18], all of which use some
source of energy, typically internal, to move or to apply forces. Such bodies,
referred to as active particles in the physics literature, mutually interact
and lead to non-equilibrium states, which may display rich spatio-temporal
evolution. The bird-flocking model of Vicsek [3], a non-equilibrium version of
a Heisenberg-spin model, is defined in discrete time for an assembly of point
particles which are distributed randomly in space. These particles try to
align with their neighbours, but with some stochastically modelled error.

Soon after the publication of the Vicsek model [3], Toner and Tu (TT)
introduced a hydrodynamic stochastic partial differential equation (PDE) that
models flocking phenomena [5, 6]. The TT velocity field obeys a generalised,
compressible, stochastically forced set of Navier-Stokes equations (NSEs),
which are not Galilean invariant. Other hydrodynamic PDEs, related to the
NSEs, were then developed to study the spatio-temporal evolution of active
fluids, such as dense bacterial suspensions [5, 6, 7, 8, 9, 10, 11, 12, 14,
15, 13, 16, 17, 18], or active nematics. In one of the simplest variants,
called either the mean-bacterial-velocity or the Toner-Tu-Swift-Hohenberg
(TTSH) model, a term consisting of the sum of a negative Laplacian and a bi-
Laplacian is added to an incompressible, deterministic TT PDE (henceforth,
ITT) [21]. For recent studies of the stochastically forced and deterministic
variants of the ITT equations we refer the reader to Refs. [22, 23, 24, 25,
26].

Although these active-matter and active-fluid PDEs have been studied
intensively over the past two decades from a physical perspective, with the
results of these investigations having undergone wide experimental comparison,
detailed methods of Navier-Stokes analysis [27, 28, 29, 30, 31, 32, 33] have
not generally been applied directly to the ITT equations111An exception is the
work of Zanger, Löwen, and Saal on the regularity of solutions of the TTSH
equations [34].. However, models related to the NSEs with an absorption term
have been studied [35], including the Brinkman-Forchheimer-extended Darcy
model of porous media [36, 35, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47].
The major difference is that these models possess a nonlinearity that breaks
the Navier-Stokes invariance enjoyed by the ITT equations.

The main aim of this paper is to consider the ITT PDEs in $d$ spatial
dimensions ($d=2,3$) using the ideas developed in Refs. [48, 49, 50, 51] in
which a combination of analysis and direct numerical simulations (DNSs) on the
$d=3$ NSEs were used to match the results of the former against those of the
latter. For the NSEs in dimension $d=3$, only weak solutions (in the sense of
Leray [27]) are known to exist. To elevate these to the status of strong
solutions, a uniqueness property would need to be proved, but this proof
remains elusive [27, 28, 29, 30, 31, 32, 33]. The first task is to show that
the ITT equations in dimension $d=3$ have equivalent weak solution properties
to their Navier-Stokes counterpart. The second task is to study the $d=2$
case, keeping in mind that the NSEs are known to be regular, i.e. global
existence and uniqueness of solutions have been proved [28]. It will be shown
that the ITT equations also turn out to be regular, although the counterpart
of bounds that remain fixed for the NSEs actually grow exponentially in time.
While numerical (pseudo-spectral) methods and experiments are able to track a
solution that evolves from specified initial conditions, methods of analysis
are unable to do this; instead, in a complementary fashion, they provide us
with constraints on solutions that evolve from all smooth initial conditions.
Thus, one should not expect the estimated bounds to be saturated, as these
take into account all smooth initial conditions, however large – see [48, 49,
50, 51]. Methods of analysis also provide upper bounds on average inverse
length scales, which can be interpreted as lower bounds on the grid sizes
necessary to resolve solutions.

In §2 we define the PDEs in dimensionless form and the quantities that are
required for our analysis. In §3 the scaling properties of the Navier-Stokes
and the ITT equations are discussed and how their similarity acts as a guide
to our choice of moments of higher derivatives of the velocity field. In §4 we
discuss energy estimates. In §5 we describe the pseudo-spectral DNS that has
been used to solve the ITT equations. §6 deals solely with the $d=2$ case : in
§6.1 we present a summary of our results, the proofs of which can be found in
B, while §6.2 and §6.4 are devoted to global regularity and numerical results
respectively. Likewise in §7 we summarise our results for $d=3$, with the
proofs found in C. In §8 we discuss the significance of our results and
compare them with similar results for related PDEs.

##  2 Dimensionless equations

The standard form of the incompressible Toner-Tu (ITT) equations is given by
[5, 6, 22] :

$\left(\partial_{t}+\lambda\mbox{\boldmath$u$}\cdot\nabla\right)\mbox{\boldmath$u$}+\nabla
p=\alpha\mbox{\boldmath$u$}+\nu\Delta\mbox{\boldmath$u$}-\beta\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2}\,.$
(2.1)

The fixed parameters $\alpha,\,\beta$ are positive and the velocity field $u$
satisfies the incompressibility condition $\mbox{div}\,\mbox{\boldmath$u$}=0$.
$\beta$ has the dimension $TL^{-2}\equiv[\nu^{-1}]$, $\alpha$ is a frequency
and $\lambda$ is a dimensionless parameter. The domain is taken to be a
periodic box $[0,\,L]^{d}_{per}$. We leave remarks until §3 on the literature
involving generalizations of this system to a
$\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2\delta}$ nonlinear term.

The first step is to introduce a typical velocity field $U_{0}$ for which we
have two definitions:

$U_{0}=\sqrt{\alpha/\beta}\qquad\mbox{and}\qquad U_{0}=\nu/L\,.$ (2.2)

Then primed dimensionless variables are defined thus :

$\displaystyle\mbox{\boldmath$x$}^{\prime}=L^{-1}\mbox{\boldmath$x$}\,;\qquad
t^{\prime}=U_{0}L^{-1}t\,;\qquad\mbox{\boldmath$u$}^{\prime}=\lambda
U_{0}^{-1}\mbox{\boldmath$u$}\,;\qquad p^{\prime}=\lambda U_{0}^{-2}p\,.$
(2.3)

This transforms (2.1) into the dimensionless ITT equations (dropping the
primes) which, from now on, will be the form used in this paper :

$\displaystyle\left(\partial_{t}+\mbox{\boldmath$u$}\cdot\nabla\right)\mbox{\boldmath$u$}+\nabla
p$ $\displaystyle=$ $\displaystyle\alpha_{0}\mbox{\boldmath$u$}+{\rm
Re}_{\nu}^{-1}\Delta\mbox{\boldmath$u$}-{\rm
Re}_{\beta}\,\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2}\,,$ (2.4)

together with the incompressibility condition
$\mbox{div}\,\mbox{\boldmath$u$}=0$. These operate on the unit periodic box
$V_{d}=[0,\,1]^{d}$. The two Reynolds numbers ${\rm Re}_{\nu}$ and ${\rm
Re}_{\beta}$ are defined as follows :

$\displaystyle{\rm Re}_{\nu}$ $\displaystyle=$
$\displaystyle\frac{U_{0}L}{\nu}\,,\qquad{\rm Re}_{\beta}=\frac{\beta
U_{0}L}{\lambda^{2}}\,,$ (2.5)

taken together with the dimensionless frequency $\alpha_{0}=L\alpha
U_{0}^{-1}$. The second choice of $U_{0}$ corresponds to ${\rm Re}_{\nu}=1$.

##  3 Invariant scaling, time averages and length scales

The incompressible NSEs possess the following well-known and powerful
invariant scaling property involving an arbitrary parameter $\ell$:

$\mbox{\boldmath$x$}^{\prime}=\ell^{-1}x;\quad
t^{\prime}=\ell^{-2}t;\quad\mbox{\boldmath$u$}=\ell^{-1}\mbox{\boldmath$u$}^{\prime};$
(3.1)

which means that these equations are valid at every scale. The effect of this
invariance is to scale the norms $\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}$
defined by

$\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}=\left(\int_{V_{d}}|\nabla^{n}\mbox{\boldmath$u$}|^{2m}dV_{d}\right)^{1/2m}$
(3.2)

in the following way:

$\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}=\ell^{-1/\alpha_{m,n,d}}\|\nabla^{{}^{\prime}n}\mbox{\boldmath$u$}^{\prime}\|_{2m}\,,$
(3.3)

where $\alpha_{n,m,d}$ is defined by222The possible confusion caused by the
labelling of the dimensionless frequency $\alpha_{0}$ and the exponents
$\alpha_{n,m}$ is unfortunate, but we continue to use it to avoid the greater
confusion of changing the notation from previous papers.

$\alpha_{n,m,d}=\frac{2m}{2m(n+1)-d}\,.$ (3.4)

The $\alpha_{n,m,d}$ are a product of the invariance property (3.1). A
dimensionless version of the norms defined in (3.3) is given by

$F_{n,m,d}:=\nu^{-1}L^{1/\alpha_{n,m,d}}\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}\,.$
(3.5)

It has been shown that, for $d=2,\,3$, and for $n\geq 1$ and $1\leq
m\leq\infty$, weak solutions of the incompressible NSEs obey [32, 33]

$\left<F_{n,m,d}^{(4-d)\alpha_{n,m,d}}\right>_{T}<\infty\,.$ (3.6)

The angular brackets $\left<\cdot\right>_{T}$ represent the time average up to
a time $T$, i.e.,

$\left<\cdot\right>_{T}=\frac{1}{T}\int_{0}^{T}\cdot\,d\tau\,.$ (3.7)

We emphasize that these brackets represent a time average, not a statistical
average. When $n=0$ then $m$ is restricted by $3<m\leq\infty$. An example
familiar to the reader is the case $n=m=1$, in which case
$(4-d)\alpha_{1,1,d}=2$ with the cancellation of the factor of $4-d$. Then
(3.6) yields the333The vorticity $\omega$ and the velocity gradient tensor
$\nabla\mbox{\boldmath$u$}$ are synonymous in $L^{2}$ when
$\mbox{div}\,\mbox{\boldmath$u$}=0$, but not in $L^{p}$ for $p>2$. familiar
bound on the time-averaged energy dissipation rate

$\varepsilon=\nu
L^{-d}\left<\int_{V}|\mbox{\boldmath$\omega$}|^{2}dV_{d}\right>_{T}\leq\nu^{3}L^{-4}{\rm
Re}^{3}.$ (3.8)

With the inverse Kolmogorov length defined by
$\lambda_{k}^{-4}=\varepsilon/\nu^{3}$, we obtain the conventional bound

$L\lambda_{k}^{-1}\leq{\rm Re}^{3/4}\,.$ (3.9)

Equation (3.6) thus expresses an infinite hierarchy of such bounds and can be
looked upon as weighted space-time averages of all derivatives of the velocity
field in every $L^{2m}$-norm. There is an informal analogy with the concept of
wavelets : higher derivatives reflect the dynamics at small scales, while
increasing the value of $m$ magnifies the larger amplitudes at each scale.

In [32, 33] it has also been shown how to define a set of inverse length
scales associated with (3.6). Consider the set of $t$-dependent length-scales
$\\{\ell_{n,m,d}(t)\\}$ defined by

$\left(L\ell_{n,m,d}^{-1}\right)^{n+1}:=F_{n,m,d}\,.$ (3.10)

This definition takes into account the scaling of the domain volume $L^{d}$
which makes (3.10) at the level of $n=m=1$ and $d=3$ consistent with the
correct definition of the energy dissipation rate used to define the
Kolmogorov length. Then we easily find that for Navier-Stokes weak solutions,
when $n\geq 1$ and $1\leq m\leq\infty$ ,

$\left<L\ell_{n,m,d}^{-1}\right>_{T}\leq
c_{n,m,d}Re^{\frac{3}{(4-d)(n+1)\alpha_{n,m,d}}}\,.$ (3.11)

When $d=3$ and $n=m=1$, then the exponent in (3.11) is
${\scriptstyle\frac{3}{4}}$, which is consistent with estimates for the
Kolmogorov length [32, 33]. Also, note that $(n+1)\alpha_{n,m,d}\to 1$ as
$n,\,m\to\infty$ .

Of course, it has been known for many years that solutions of the NSEs in
$d=2$ dimensions are regular [28, 29], but expressing (3.6) in integer
$d$-dimensions $d=1,\,2,\,3$ rolls together into one line all the known two-
and three-dimensional Navier-Stokes solution results, such as the class of
weak solution $d=3$ time averages found by Foias, Guillopé and Temam [31] in
their pioneering paper in 1981. It has been explained in Ref. [32] that, for a
full existence and uniqueness proof in the $d=3$ case, a factor of
$2\alpha_{n,m,3}$ would be needed in the exponent in (3.6). However, no
evidence exists for the existence of bounds with this necessary factor of $2$.
It is possible that the Leray-Hopf weak solutions are all that exist.

By inspection it is clear that the ITT equations respect the invariant scaling
possessed by the NSEs equations, apart from the linear-pumping term. However,
there is a significant literature on a more general class of equations where
the $\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2}$-term is replaced by
$\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2\delta}$, which is the case in the
Brinkman-Forchheimer extended Darcy model arising in porous media. The paper
by Titi and Trabelsi [36] contains a wide literature survey, but we also refer
the reader to [35, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]. When
$\delta>1$ , the invariant scaling property of the NSEs is broken. This leads
to the bounding of time-averaged norms, higher than those available to ITT,
which eventually lead to the regularity of solutions in the $d=3$ case. The
work in our paper is different in two important respects. Firstly, we focus on
the critical value $\delta=1$ and thus remain faithful to the scaling property
in (3.1). Secondly, the ITT equations in (2.4) have a linear term
$\alpha_{0}\mbox{\boldmath$u$}$ which, while trivial in a purely functional
setting, is nevertheless physically important in the creation of equilibrated
or frozen states, which appear to dominate the dynamics in our DNSs.

The parallel scaling properties of the ITT equations and the NSEs suggest that
the exponents $\alpha_{n,m,d}$ in (3.4) should be the same in both cases.
Therefore, taking into account the factor of $4-d$ in the exponent in the
$d=2$ case, we define

$P_{n,m}=\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}^{2\alpha_{n,m,2}}\,,\qquad(d=2)\,,$
(3.12)

where

$\alpha_{n,m,2}=\frac{m}{m(n+1)-1}$ (3.13)

and a set of inverse length scales equivalent to (3.11)

$\left<L\ell_{n,m,2}^{-1}\right>_{T}\leq
c_{n,m,2}\left<P_{n,m}\right>_{T}^{\frac{1}{2(n+1)\alpha_{n,m,2}}}\,.$ (3.14)

When $d=3$

$Q_{n,m}=\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}^{\alpha_{n,m,3}}\,,\qquad(d=3)\,,$
(3.15)

where

$\alpha_{n,m,3}=\frac{2m}{2m(n+1)-3}$ (3.16)

and a set of inverse length scales equivalent to (3.11)

$\left<L\ell_{n,m,3}^{-1}\right>_{T}\leq
c_{n,m,3}\left<Q_{n,m}\right>_{T}^{\frac{1}{(n+1)\alpha_{n,m,3}}}\,.$ (3.17)

General bounds on $\left<P_{n,m}\right>_{T}$, expressed as a function of
$\alpha_{0}$, ${\rm Re}_{\nu}$ and $\mathcal{A}_{0}$, are estimated in §6. For
$d=3$, a narrower class of bounds on $\left<Q_{n,m}\right>_{T}$ has been given
in §7.

##  4 Energy estimates

In keeping with standard Navier-Stokes notation, we define $n$ derivatives of
$u$ in $L^{2}(V_{d})$ as

$H_{n}=\int_{V_{d}}|\nabla^{n}\mbox{\boldmath$u$}|^{2}dV_{d}\,.$ (4.1)

Given the close relationship between the ITT equations and the incompressible
NSEs, a formal approach is taken following the method used by [31]. In this
notation the energy $H_{0}$ and and enstrophy $H_{1}$ are

$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle\int_{V_{d}}|\mbox{\boldmath$u$}|^{2}dV_{d}\,;$ (4.2)
$\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle\int_{V_{d}}|\nabla\mbox{\boldmath$u$}|^{2}dV_{d}=\int_{V_{d}}|\mbox{\boldmath$\omega$}|^{2}dV_{d}\,.$
(4.3)

Using standard methods [27, 28], a Leray-type energy inequality is easily
derived

${\scriptstyle\frac{1}{2}}\dot{H}_{0}+{\rm Re}_{\nu}^{-1}H_{1}+{\rm
Re}_{\beta}\int_{V}|\mbox{\boldmath$u$}|^{4}dV_{d}\leq\alpha_{0}H_{0}$ (4.4)

from which we drop the $H_{1}$-term444Poincaré’s inequality cannot be applied
because, unlike the NSEs, the spatial average of $u$ is not zero. and apply a
Hölder inequality to the $L^{4}$-term to produce a simple differential
inequality for $H_{0}$

${\scriptstyle\frac{1}{2}}\dot{H}_{0}\leq\alpha_{0}H_{0}-{\rm
Re}_{\beta}H_{0}^{2}\,.$ (4.5)

Thus, equilibration of the right hand side occurs at

$H_{0,\,equil}=\alpha_{0}{\rm Re}_{\beta}^{-1}\,:=\,\mathcal{A}_{0}\,,$ (4.6)

where we designate $\mathcal{A}_{0}$ as the activity parameter. By using the
time-average definition in (3.7), from (4.4) we find the following estimates

$\displaystyle\left<H_{0}\right>_{T}$ $\displaystyle\leq$
$\displaystyle\mathcal{A}_{0}\,,$ (4.7) $\displaystyle\left<H_{1}\right>_{T}$
$\displaystyle\leq$ $\displaystyle\alpha_{0}\mathcal{A}_{0}{\rm Re}_{\nu}\,,$
(4.8)

together with the time average of both the $L^{4}$-norm and the ratio of
$H_{1}$ to $H_{0}$

$\left<\int_{V}|\mbox{\boldmath$u$}|^{4}dV_{d}\right>_{T}\leq\mathcal{A}_{0}^{2}\,,\qquad\qquad\left<\frac{H_{1}}{H_{0}}\right>_{T}\leq\alpha_{0}{\rm
Re}_{\nu}\,.$ (4.9)

The inequalities (4.7)–(4.9) each have an $O\left(T^{-1}\right)$ correction
term that will be dropped from now on. These results are true in every
dimension. With respect to initial data on $H_{0}$ :

  1. 1.

For initial data $H_{0}(0)>\mathcal{A}_{0}$ , the sign of $\dot{H}_{0}$ in
this region is negative and thus $H_{0}$ decreases to $\mathcal{A}_{0}$.

  2. 2.

For initial data $H_{0}(0)<\mathcal{A}_{0}$ , $\dot{H}_{0}$ can take either
sign. If $\dot{H}_{0}(t)>0$ then $H_{0}$ will grow to reach $\mathcal{A}_{0}$
(a frozen state) but cannot pass through it. If $\dot{H}_{0}(t)<0$ then
$H_{0}$ will decay.

The choice of $n=m=1$ makes $(4-d)\alpha_{1,1,d}=2$, whereupon the factor of
$4-d$ cancels. Thus

$\left<P_{1,1}\right>_{T}\leq\alpha_{0}\mathcal{A}_{0}{\rm Re}_{\nu}\,,$
(4.10)

and

$\left<Q_{1,1}\right>_{T}\leq\alpha_{0}\mathcal{A}_{0}{\rm Re}_{\nu}\,.$
(4.11)

##  5 Numerical Methods

Run | $d$ | $N$ | $\delta t$ | $\nu$ | $\alpha$ | $\beta$
---|---|---|---|---|---|---
A1 | 2 | 2048 | $2\times 10^{-4}$ | $2.87\times 10^{-1}$ | $10\times 10^{1}$ | $5$
A2 | 2 | 2048 | $2\times 10^{-4}$ | $1.41\times 10^{-1}$ | $10\times 10^{1}$ | $5$
A3 | 2 | 2048 | $2\times 10^{-4}$ | $7.07\times 10^{-2}$ | $10\times 10^{1}$ | $5$
A4 | 2 | 2048 | $2\times 10^{-4}$ | $3.53\times 10^{-2}$ | $10\times 10^{1}$ | $5$
A5 | 2 | 2048 | $2\times 10^{-4}$ | $2.36\times 10^{-2}$ | $10\times 10^{1}$ | $5$
A6 | 2 | 2048 | $2\times 10^{-4}$ | $1.77\times 10^{-2}$ | $10\times 10^{1}$ | $5$
A7 | 2 | 2048 | $2\times 10^{-4}$ | $1\times 10^{-2}$ | $10\times 10^{1}$ | $5$
A8 | 2 | 2048 | $2\times 10^{-4}$ | $8.8\times 10^{-2}$ | $10\times 10^{1}$ | $5$
F1 | 2 | 2048 | $2\times 10^{-4}$ | $6.2\times 10^{-1}$ | $1$ | $1$
F2 | 2 | 2048 | $2\times 10^{-4}$ | $1.2\times 10^{-1}$ | $1$ | $1$
F3 | 2 | 2048 | $2\times 10^{-4}$ | $6.0\times 10^{-2}$ | $1$ | $1$
F4 | 2 | 2048 | $2\times 10^{-4}$ | $3.0\times 10^{-2}$ | $1$ | $1$
F5 | 2 | 2048 | $2\times 10^{-4}$ | $1.5\times 10^{-2}$ | $1$ | $1$
F6 | 2 | 2048 | $2\times 10^{-4}$ | $7.0\times 10^{-3}$ | $1$ | $1$
F7 | 2 | 2048 | $2\times 10^{-4}$ | $3.1\times 10^{-3}$ | $1$ | $1$
B1 | 3 | 512 | $1\times 10^{-3}$ | $5\times 10^{-1}$ | $1\times 10^{1}$ | $1\times 10^{-1}$
B2 | 3 | 512 | $1\times 10^{-3}$ | $5\times 10^{-2}$ | $1\times 10^{1}$ | $1\times 10^{-1}$
B3 | 3 | 512 | $5\times 10^{-4}$ | $1\times 10^{-2}$ | $1\times 10^{1}$ | $1\times 10^{-1}$
Table 1:  The parameters for our DNSs: $d$ is the dimension, $N^{d}$ the
number of collocation points, and $\delta$ the time step.For all our runs,
$\lambda=1$. Given these parameters, the Reynolds numbers follow from Eqs.
(2.2) and (2.5). Parameters for other runs are given in the Supplemental
Material.

For our DNS of the $d$-dimensional ITT Eq. (2.1), we use a Fourier
pseudospectral method [54] on periodic domains (a square in $d=2$ and a cube
in $d=3$), with sides of length $L=2\pi$, and $N^{d}$ collocation points. We
employ the second-order exponential time-difference scheme, ETDRK2, for time
evolution in Fourier-space [55]. We list the parameters for our DNS runs in
Table 1. Our simulations are numerically robust insofar as the Courant-
Friedrichs-Lewy number, $C<C_{max}$, where
$C=\left(\sum_{i=1}^{i=d}U_{i}/h_{i}\right)\delta t$ with
$U_{i}=~{}\sup_{x}u_{i}(\mbox{\boldmath$x$},t)$ and $h_{i}$ the minimum grid
spacing between two collocation points. Typically $C_{max}=1$ for our explicit
numerical schemes and for all our simulations (Table 1) $C\simeq 0.62\pm
0.04$.

The dimensional version of the ITT equations (2.1) has four parameters
$\lambda,\,\alpha,\,\nu$ and $\beta$ which reduce to the three dimensionless
numbers $Re_{\nu},\,Re_{\beta}$ and $\alpha_{0}$ in the non-dimensionalized
version (2.4). As explained in (2.2), we have found it convenient to define
the typical velocity field $U_{0}$ in two particular ways :
$U_{0}=\sqrt{\alpha/\beta}$ and $U_{0}=\nu/L$. The latter case restricts
$Re_{\nu}$ to the value $Re_{\nu}=1$ but allows us to explore a more diverse
range of $\alpha_{0}$ and $Re_{\beta}$.

##  6 Summary of results in the $d=2$ case

The methods used in the analysis sections of this paper are based on the
differential inequalities explained in A. The proof of the results in the
following subsections are given in B. Within these estimates, various
multiplicative constants $c,\,c_{m}$ and $c_{n,m}$ appear, which should be
read as generic constants that may differ from line to line. These constants
are algebraic in $n,m$ but are not usually given explicitly : see A. We remark
that none of the bounds displayed in the following sections are saturated,
although more drastic initial conditions might get closer to this state. For
simplicity, it is also assumed that ${\rm Re}_{\nu}\gg{\rm Re}_{\beta}$. In
the choice $U_{0}=\nu/L$ where ${\rm Re}_{\nu}=1$, in which case estimates can
be re-calculated from the material in the Appendices.

###  6.1 Estimates for $\left<P_{n,m}\right>_{T}$

By using the definition of $P_{n,m}$ in (3.12) as our guide, for $m=1$, we
have $\alpha_{n,1,2}=1/n$ , so

$\left<P_{n,1}\right>_{T}=\left<H_{n}^{1/n}\right>_{T}\,.$ (6.1)

The estimate for $\left<P_{1,1}\right>_{T}$ in (4.10) can be used to compute a
series of other inequalities. The proofs can be found in B. Below is a summary
:

(i) Firstly, we wish to estimate $\left<P_{n,1}\right>_{T}$ for $n\geq 1$. In
inequality (B.8) in B.1, for $n=2$ it is shown that

$\left<P_{2,1}\right>_{T}\leq c\,\alpha_{0}\left(\alpha_{0}\mathcal{A}_{0}{\rm
Re}_{\nu}^{3}\right)^{1/2}\,.$ (6.2)

(ii) More generally, in inequality (B.22) B.4, it is shown that for $n\geq 2$

$\displaystyle\left<P_{n,1}\right>_{T}$ $\displaystyle\leq$ $\displaystyle
c_{n,1}\alpha_{0}^{2/n}\left(\alpha_{0}\mathcal{A}_{0}{\rm
Re}_{\nu}^{3}\right)^{\frac{n-1}{n}}\,.$ (6.3)

(iii) In inequality (B.11) in B.2 it is shown that

$\left<P_{1,m}\right>_{T}=\left<\|\nabla\mbox{\boldmath$u$}\|_{2m}^{\frac{2m}{2m-1}}\right>_{T}\,,$
(6.4)

displayed in row 3 of Fig. 4, satisfies

$\displaystyle\left<P_{1,m}\right>_{T}$ $\displaystyle\leq$ $\displaystyle
c_{m}\left(\alpha_{0}{\rm
Re}_{\nu}\right)^{\frac{3m-2}{2m-1}}\mathcal{A}_{0}^{\frac{m}{2m-1}}\,.$ (6.5)

(B.12) shows that in the limit $m\to\infty$, we also have

$\left<\|\nabla\mbox{\boldmath$u$}\|_{\infty}\right>_{T}\leq
c\,\left(\alpha_{0}{\rm Re}_{\nu}\right)^{3/2}\mathcal{A}_{0}^{1/2}\,.$ (6.6)

(iv) With the definition

$\left<P_{0,m}\right>_{T}=\left<\|\mbox{\boldmath$u$}\|_{2m}^{\frac{2m}{m-1}}\right>_{T}\,,$
(6.7)

from (B.15) in B.3, it is shown that, for $m>2$ ,

$\left<P_{0,m}\right>_{T}\leq
c\,\mathcal{A}_{0}^{\frac{m}{m-1}}\left(\alpha_{0}{\rm
Re}_{\nu}\right)^{\frac{m-2}{m-1}}\,.$ (6.8)

(B.16) shows that, in the limit $m\to\infty$, we have

$\left<P_{0,\infty}\right>_{T}=\left<\|\mbox{\boldmath$u$}\|_{\infty}^{2}\right>_{T}\leq
c\,\alpha_{0}\mathcal{A}_{0}{\rm Re}_{\nu}\,.$ (6.9)

(v) In inequality (B.29) in B.5 it is shown that, for $n\geq 2$,

$\left<P_{n,m}\right>_{T}\leq
c_{n,m}\alpha_{0}^{\frac{2m}{m(n+1)-1}}\left(\alpha_{0}\mathcal{A}_{0}{\rm
Re}_{\nu}^{3}\right)^{\frac{mn-1}{m(n+1)-1}}\,.$ (6.10)

###  6.2 Regularity : exponential bounds in $d=2$ dimensions

The difference between the $2d$ and $3d$ NSEs lies in the absence of the
vortex stretching term
$\mbox{\boldmath$\omega$}\cdot\nabla\mbox{\boldmath$u$}$, which is zero when
$d=2$. For the ITT equations (2.4), the evolution equation for the vorticity
is

$\left(\partial_{t}+\mbox{\boldmath$u$}\cdot\nabla\right)\mbox{\boldmath$\omega$}=\alpha_{0}\mbox{\boldmath$\omega$}+{\rm
Re}_{\nu}^{-1}\Delta\mbox{\boldmath$\omega$}-{\rm
Re}_{\beta}\,\mbox{curl}\left(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2}\right)$
(6.11)

in which the
$\mbox{curl}\left(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2}\right)$-term
appears to create another form of vortex stretching. Specifically we have

$\displaystyle{\scriptstyle\frac{1}{2}}\dot{H}_{1}$ $\displaystyle=$
$\displaystyle\alpha_{0}H_{1}-{\rm Re}_{\nu}^{-1}H_{2}-{\rm
Re}_{\beta}\int_{V}\mbox{\boldmath$\omega$}\cdot\mbox{curl}(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2})dV$
(6.12) $\displaystyle\leq$
$\displaystyle\alpha_{0}H_{1}-{\scriptstyle\frac{1}{2}}{\rm
Re}_{\nu}^{-1}H_{2}+{\scriptstyle\frac{1}{2}}{\rm Re}_{\beta}^{2}{\rm
Re}_{\nu}\int_{V}|\mbox{\boldmath$u$}|^{6}dV\,,$

where we have integrated by parts and have then used a Hölder inequality. Then
we use a Gagliardo-Nirenberg inequality

$\|\mbox{\boldmath$u$}\|_{6}\leq
c\,\|\nabla\mbox{\boldmath$u$}\|_{2}^{{\scriptstyle\frac{1}{3}}}\|\mbox{\boldmath$u$}\|_{4}^{{\scriptstyle\frac{2}{3}}}\,,$
(6.13)

which can be re-expressed as

$\int_{V}|\mbox{\boldmath$u$}|^{6}dV\leq
c\,H_{1}\int_{V}|\mbox{\boldmath$u$}|^{4}dV\,.$ (6.14)

Thus we can write

${\scriptstyle\frac{1}{2}}\dot{H}_{1}\leq\left(\alpha_{0}+{\scriptstyle\frac{1}{2}}{\rm
Re}_{\beta}^{2}{\rm
Re}_{\nu}\int_{V}|\mbox{\boldmath$u$}|^{4}\,dV\right)H_{1}\,,$ (6.15)

whence we can use the extra piece of information afforded to us by the bounded
time integral expressed in (4.9). Thus (6.15) can be integrated to become

$\displaystyle H_{1}(T)$ $\displaystyle\leq$ $\displaystyle
H_{1}(0)\exp\left\\{\int_{0}^{T}\left(\alpha_{0}+c\,{\rm Re}_{\beta}^{2}{\rm
Re}_{\nu}\|\mbox{\boldmath$u$}\|_{4}^{4}\right)\,d\tau\right\\}$ (6.16)
$\displaystyle\leq$ $\displaystyle
H_{1}(0)\exp\left\\{\alpha_{0}\left(1+c\,{\rm Re}_{\beta}^{2}{\rm
Re}_{\nu}\mathcal{A}_{0}^{2}\right)T\right\\}\,,$

which is finite for every finite $T$. Control over the $H_{1}$-norm
establishes global regularity in this $2d$ case but not a global attractor,
which requires a uniform bound for all $t$.

###  6.3 Spectral energy budget

We discuss the role of contributions to the energy from the nonlinear terms in
the ITT equations. The shell-averaged energy spectrum is defined by

$\mathcal{E}(k)=\frac{1}{2}\sum_{k^{\prime}=k-1/2}^{k^{\prime}=k+1/2}\left<\widetilde{\mbox{\boldmath$u$}}(\mbox{\boldmath$k$}^{\prime})\cdot\widetilde{\mbox{\boldmath$u$}}(-\mbox{\boldmath$k$}^{\prime})\right>_{t}\,,$
(6.17)

where a tilde denotes a spatial Fourier transform and the wave vectors $k$ and
$\mbox{\boldmath$k$}^{\prime}$ have moduli $k$ and $k^{\prime}$, respectively.
In Figs. 1 and 2 we give illustrative plots of the time series of the total
kinetic energy

$E_{tot}(t)\equiv{\scriptstyle\frac{1}{2}}L^{-d}\int_{V_{d}}|\mbox{\boldmath$u$}|^{2}\,dV_{d}$
(6.18)

together with $\mathcal{E}(k)$, filled contour plots of the vorticity (in
$d=2$) and isosurfaces of the modulus of the vorticity (in $d=3$). From Fig.
1(a) we surmise the existence of temporally frozen states in $d=2$, because
$dE_{tot}/dt\simeq 0$ after the period of initial transients. Both frozen and
turbulent states can also be characterized by the filled contour plots of the
vorticity presented in column 3 of Fig. 1; these states depend on $\nu$,
$\alpha$, and $\beta$, but not on the initial conditions ; a change in the
initial conditions changes the time for which initial transients last. We also
find statistically steady states, with fluctuations in the total energy, in
both $d=2$ (Fig. 1(d)) and $d=3$ (Fig. 2(a)), similar to those observed in
homogeneous isotropic fluid turbulence.

The $k$-dependent energy budget is given by

${}\partial_{t}\mathcal{E}(k)=T(k)-T_{\beta}(k)+T_{\alpha}(k)-2\Gamma_{0}k^{2}\mathcal{E}(k)$
(6.19)

where $T_{\alpha}(k)$, the spectral energy contribution from the $\alpha$
term, and $T(k)$ and $T_{\beta}(k)$ are, respectively, the $k$-shell averaged
contributions from the advective and cubic terms in the ITT equations, are:

$\displaystyle{}T(k)$
$\displaystyle=-\lambda\sum_{k^{\prime}=k-1/2}^{k^{\prime}=k+1/2}\left<\widetilde{\mbox{\boldmath$u$}}(-\mbox{\boldmath$k$})\cdot\textbf{P}\left(\mbox{\boldmath$k$}^{\prime}\right)\cdot\widetilde{\left(\mbox{\boldmath$u$}\cdot\nabla\mbox{\boldmath$u$}\right)}(\mbox{\boldmath$k$}^{\prime})\right>_{t}\,;$
(6.20) $\displaystyle T_{\beta}(k)$
$\displaystyle=\beta\sum_{k^{\prime}=k-1/2}^{k^{\prime}=k+1/2}\left<\widetilde{\mbox{\boldmath$u$}}(-\mbox{\boldmath$k$})\cdot\textbf{P}(\mbox{\boldmath$k$}^{\prime})\cdot\widetilde{(|\mbox{\boldmath$u$}|^{2}\mbox{\boldmath$u$})}(\mbox{\boldmath$k$}^{\prime})\right>_{t}\,;$
(6.21) $\displaystyle T_{\alpha}$ $\displaystyle=-2\alpha\mathcal{E}(k)\,;$
(6.22)

here,

$P_{ij}(\mbox{\boldmath$k$})=\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}$ (6.23)

is the transverse projector and $\left<\cdot\right>_{t}$ the average over time
$t$. The flux of energy arising from the advective and the cubic terms are,
respectively:

$\displaystyle{}\Pi(k)$
$\displaystyle=-\sum_{k^{\prime}=0}^{k^{\prime}=k}T(k^{\prime})\,;$ (6.24)
$\displaystyle\Pi_{\beta}(k)$
$\displaystyle=-\sum_{k^{\prime}=0}^{k^{\prime}=k}T_{\beta}(k^{\prime})\,.$
(6.25)

In Fig. 3(a), we plot versus $k$, $T_{\beta}(k)$ and $T_{\alpha}(k)$; we omit
the points $T_{\alpha}(1)$ and $T_{\beta}(1)$ as
$T_{\alpha}(1)>>T_{\alpha}(k)$ and $T_{\beta}(1)<<T_{\beta}(k)$. We note from
$T_{\beta}(k)$ that the cubic term is an energy sink and is dominant at
small-$k$ modes; by contrast, the $\alpha$ term acts as an energy source.

In Figs. 3(b) and (c) we give plots versus $k$ of the fluxes $\Pi(k)$ and
$\Pi_{\beta}(k)$, respectively. The flux associated with the cubic term
$\Pi_{\beta}(k)<0$, as we expect from an energy sink. However, the advective
term is neither a sink nor a source of energy. Therefore, the total area under
the curve in Fig. 3 (a) is zero. We also note that there is no region of $k$
for which $\Pi(k)\simeq$ a constant. This indicates the absence of a
conventional inertial range.

(a)(b)(c)(d)(e)(f)

Figure 1: Plots for runs F7 (row 1) and A6 (row2): column (1) contains plots
versus the time $t$ of the total energy $E_{tot}(t)$; column (2) contains log-
log plots versus $k$ of the energy spectrum $E(k)$; the filled contour plots
in column (3) are of $\omega$ at a representative time.

(a)(b)(c)

Figure 2: (Colour online) For run B2; plot versus $t$ of the total kinetic
energy $E_{tot}(t)$ (column 1, row 1); log-log plot versus $k$ of the energy
spectrum $E(k)$ (column 1, row2); iso-surfaces of the modulus of the vorticity
field (column 2).

(a)(b)(c)

Figure 3: (Color online) Semi-log plots versus $k$ of (a) $T_{\alpha}(k)$
(blue) and $T_{\beta}(k)$ (red), (b) $\Pi(k)$ (red), and (c) $\Pi_{\beta}(k)$
(blue), for an illustrative simulation (A4) from Table. 1

###  6.4 Numerical results for $P_{n,m}$

We now present plots in Fig. 4 of $P_{n,m}(t)$ versus time $t$ for the two
values $n=0$ and $n=1$, with a sequence of values of $m=1,\,\ldots,\,10$.
Hölder’s inequality insists that, for fixed $n$, the norms $\|\cdot\|_{2m}$
must be ordered with increasing $m$, such that
$\|\cdot\|_{2m}\leq\|\cdot\|_{2(m+1)}$; but the $\alpha_{n,m,d=2}$ decrease as
$m$ increases. Thus, it is technically possible for the $P_{n,m}$ to be
ordered either way : i.e., an increasing regime $P_{n,m}\leq P_{n,m+1}$ or a
decreasing regime $P_{n,m}\geq P_{n,m+1}$. The latter regime was originally
observed numerically for the $d=3$ NSEs [48, 49] and was also discussed in
[50], although no obvious reason for this particular ordering was deduced.
Moreover, no crossing of curves that represented different values of $m$ has
been observed.

For the $d=2$ case, when $U_{0}=\sqrt{\alpha/\beta}$ (panel A), we have
observed both regimes, but only the decreasing regime $P_{n,m}\geq P_{n,m+1}$
when $U_{0}=\nu/L$ (panel B). This is illustrated in Fig. 4 for run A6 : In
panel A ($U_{0}=\sqrt{\alpha/\beta}$) the plots of $P_{1,m}$ can cross each
other at different times, as we can see clearly in the expanded plots in the
second rows ; such crossings do not occur in panel B ($U_{0}=\nu/L$).
Furthermore, the plots versus $m$ of $\left<P_{1,m}\right>_{T}$ (third row)
decreases monotonically with increasing $m$ in panel B but not in panel A.

In Fig. 5 we present plots for $U_{0}=\sqrt{\alpha/\beta}$ (panel A) and
$U_{0}=\nu/L$ (panel B), runs A1-A8, to illustrate whether the bound in 4.10
is saturated : In the first row we plot $\left<P_{1,1}\right>_{T}$ (solid
black line) versus $Re_{\nu}$ (panel A) and $\alpha_{0}$ (panel B) ; the black
dashed line denotes $Re_{\nu}\,\alpha_{0}\,\mathcal{A}_{0}$, which is the
right-hand side (RHS) of (4.10) . In the second row we present plots versus
$Re_{\nu}$ (panel A) and $\alpha_{0}$ (panel B) of $\left<P_{0,m}\right>_{T}$
and $\left<P_{1,m}\right>_{T}$, for $m=2,\,\ldots\,,10$. Note that curves for
$\left<P_{1,m}\right>_{T}$ can cross as $Re_{\nu}$ increases (panel A) ; by
contrast, they do not cross as $\alpha_{0}$ increases (panel B). Similar plots
for other representative runs are given in the Supplemental Material.

AB Figure 4: (colour online) Illustrative plots for
$U_{0}=\sqrt{\alpha/\beta}$ (panel A) and $U_{0}=\nu/L$ (panel B) for run A6
(see Table 1): First and second rows : plots versus $t$ of $P_{0,m}$ and
$P_{1,m}$ ; the plots in the second row are expanded versions of small
segments of the plots in the first row. Third row : Plots versus $m$ of
$\left<P_{0,m}\right>_{T}$ and $\left<P_{1,m}\right>_{T}$. Curves for
$m=2,3,4,5,6,7,8,9$, and $10$ are drawn in red, pink; violet, green, cyan,
maroon, blue, orange, and yellow, respectively. Similar plots for other
representative runs are given in the Supplemental Material. AB Figure 5:
(Colour online) Illustrative plots for $U_{0}=\sqrt{\alpha/\beta}$ (panel A)
and $U_{0}=\nu/L$ (panel B) for $d=2$, runs A1-A8 (see Table 1) : First row :
plots versus $Re_{\nu}$ (panel A) and $\alpha_{0}$ (panel B) of
$\left<P_{1,1}\right>_{T}$ ( solid black line) and
$Re_{\nu}\,\alpha_{0}\,\mathcal{A}_{0}$ (dashed black line). Second row :
Plots versus $Re_{\nu}$ (panel A) and $\alpha_{0}$ (panel B) of
$\left<P_{0,m}\right>_{T}$ and $\left<P_{1,m}\right>_{T}$. Curves for
$m=2,3,4,5,6,7,8,9$, and $10$ are drawn in red, pink; violet, green, cyan,
maroon, blue, orange, and yellow, respectively. Similar plots for other
representative runs are given in the Supplemental Material.

##  7 Summary of results in the $d=3$ case

The proof of the results in the following subsections are given in C. The
methods used there are based on the differential inequalities explained in A.

###  7.1 Estimates for $\left<Q_{n,m}\right>_{T}$

Results in the $d=3$ case are more restricted, which reflects the open status
of the regularity problem. Nevertheless, time averages of various $Q_{n,m}$ of
Navier-Stokes type can be found [32, 33]. In addition to a bound on
$\left<Q_{1,1}\right>_{T}$ , as in (4.11), our results from C are summarised
thus : from (C.16) we have

$\left<Q_{2,1}\right>_{T}\leq c\,\alpha_{0}{\rm Re}_{\nu}^{2}.$ (7.1)

We also find that for $n\geq 2$ and $m\geq 1$,

$\left<Q_{n,m}\right>_{T}<\infty\,$ (7.2)

although estimating the right hand side is a difficult calculation that we
have omitted (see (C.23)).

Moreover, with

$Q_{0,m}=\|\mbox{\boldmath$u$}\|_{2m}^{\frac{2m}{2m-3}}\,,$ (7.3)

for $m>2$, (C.28) shows that

$\left<Q_{0,m}\right>_{T}\leq
c\,\mathcal{A}_{0}^{\frac{2(m+3)}{5(2m-3)}}\left(\alpha_{0}{\rm
Re}_{\nu}^{2}\right)^{\frac{9(m-2)}{5(2m-3)}}\,.$ (7.4)

(C.29) also shows that, in the limit $m\to\infty$,

$\displaystyle\left<\|\mbox{\boldmath$u$}\|_{\infty}\right>_{T}$
$\displaystyle\leq$ $\displaystyle
c\,\alpha_{0}^{9/10}\mathcal{A}_{0}^{1/5}{\rm Re}_{\nu}^{9/5}\,.$ (7.5)

###  7.2 Numerical results for $Q_{n,m}$

AB Figure 6: (Colour online) Illustrative plots for
$U_{0}=\sqrt{\alpha/\beta}$ (panel A) and $U_{0}=\nu/L$ (panel B) for $d=3$,
run B2 (see Table 1): First and second rows: plots versus $t$ of $Q_{0,m}$ and
$Q_{1,m}$; the plots in the second row are expanded versions of small segments
of the plots in the first row. Third row: Plots versus $m$ of
$\left<Q_{0,m}\right>_{T}$ and $\left<Q_{1,m}\right>_{T}$. Curves for
$m=2,3,4,5,6,7,8,9$, and $10$ are drawn in red, pink; violet, green, cyan,
maroon, blue, orange, and yellow, respectively. Similar plots for other
representative runs are given in the Supplemental Material. AB Figure 7:
(Colour online) Illustrative plots for $U_{0}=\sqrt{\alpha/\beta}$ (panel A)
and $U_{0}=\nu/L$ (panel B) for $d=3$, runs A1-A8 (see Table 1): First row :
plots versus $Re_{\nu}$ (panel A) and $\alpha_{0}$ (panel B) of
$\left<Q_{1,1}\right>_{T}$ ( solid black line) and
$Re_{\nu}\alpha_{0}\mathcal{A}_{0}$ (dashed black line). Second row : Plots
versus $Re_{\nu}$ (panel A) and $\alpha_{0}$ (panel B) of
$\left<Q_{0,m}\right>_{T}$ and $\left<Q_{1,m}\right>_{T}$. Curves for
$m=1,2,3,4,5,6,7,8,9$, and $10$ are drawn in black, red, pink, violet, green,
cyan, maroon, blue, orange, and yellow, respectively. Similar plots for other
representative runs are given in the Supplemental Material.

Fig. 6 shows plots of $Q_{n,m}(t)$ versus time $t$ for the two values $n=0$
and $n=1$, with a sequence of values of $m=1,\,\ldots,\,10$. Again, as in the
$d=2$ case, there are two regimes, namely, an increasing regime $Q_{n,m}\leq
Q_{n,m+1}$ and a decreasing regime $Q_{n,m}\geq Q_{n,m+1}$, because the norms
$\|\cdot\|_{2m}$ must be ordered with increasing $m$, such that
$\|\cdot\|_{2m}\leq\|\cdot\|_{2(m+1)}$; but the $\alpha_{n,m,d=3}$ decrease as
$m$ increases.

When $U_{0}=\sqrt{\alpha/\beta}$ (panel A) we have observed both these
regimes, but only the decreasing regime $Q_{n,m}\geq Q_{n,m+1}$ when
$U_{0}=\nu/L$ (panel B) has been used. We illustrate this with plots in Fig. 6
for run B2 : In panel A ($U_{0}=\sqrt{\alpha/\beta}$) the plots of $Q_{1,m}$
can cross each other at different times, as we can see clearly in the expanded
plots in the second row; such crossings do not occur in panel B
($U_{0}=\nu/L$). Furthermore, the plots versus $m$ of
$\left<Q_{1,m}\right>_{T}$ (third row) decreases monotonically with increasing
$m$ in panel B but not in panel A. Note, in particular, that $Q_{1,m}$ is
almost equivalent to a non-dimensionalized version of the
$D_{m}=\|\mbox{\boldmath$\omega$}\|^{\alpha_{n,m,d=3}}_{2m}$ introduced in
[48] and plotted there for the $d=3$ NSEs (see the Supplemental Material). The
only difference here is that we are plotting
$\|\nabla\mbox{\boldmath$u$}\|_{2m}$ and not
$\|\mbox{\boldmath$\omega$}\|_{2m}$ : the two are identical only when $m=1$.

In Fig. 7 we present plots for $U_{0}=\sqrt{\alpha/\beta}$ (panel A) and
$U_{0}=\nu/L$ (panel B), runs B1-B3, to illustrate whether the bound in (4.10)
is saturated : In the first row we plot $\left<Q_{1,1}\right>_{T}$ (panel B) ;
the black dashed line denotes $Re_{\nu}\,\alpha_{0}\,\mathcal{A}_{0}$, which
is the RHS of (4.10). In the second row we present plots versus $Re_{\nu}$
(panel A) and $\alpha_{0}$ (panel B) of $\left<Q_{0,m}\right>_{T}$ and
$\left<Q_{1,m}\right>_{T}$, for $m=2,\,\ldots\,,10$. Note that curves for
$\left<Q_{1,m}\right>_{T}$ can cross as $Re_{\nu}$ increases (panel A) ; by
contrast, they do not cross as $\alpha_{0}$ increases (panel B). Similar plots
for other representative runs are given in the Supplemental Material.

##  8 Conclusions

We have married the two approaches of the analysis of solutions of the ITT
equations through the estimation of weighted time averages, together with the
results of numerical simulations. To achieve this we have invoked the similar
scaling properties between the ITT equations and the NSEs : see §3. There are,
however, two important differences. Usually the NSEs are considered either in
the decaying or the additively forced case, whereas the ITT equations have no
additive forcing but instead possess a linear-activity term
$\alpha_{0}\mbox{\boldmath$u$}$ which, in effect, pumps energy into the
system. Dynamically the effect of this term, together with the negative cubic
term, creates a platform for either temporally frozen solutions or
statistically steady states. Furthermore, it is shown in §6.2 that the
greatest contrast with the $2d$ NSEs lies at the level of the vorticity
equation. Whereas the absence of vortex stretching in $2d$ NSEs lies at the
root of the regularity of its solutions, the
$\mbox{curl}\left(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2}\right)$-term in
the ITT equations appears to recreate another form of vortex stretching.
However, the extra piece of information afforded to us is the bounded time
average of the $L^{4}$ norm of the velocity field expressed in (4.9), which is
just enough to recover regularity, but only to the degree that bounds are
exponential in time – see (6.11). Thus, we fall just short of proving the
existence of a global attractor. Results that lie in parallel with those of
the NSEs in both spatial dimensions are the existence of bounded infinite
hierarchies of time averages : i.e., estimates of $\left<P_{n,m}\right>_{T}$
and $\left<Q_{n,m}\right>_{T}$ , whose bounds are calculated in B and C and
summarised in §6 and §7, together with the results of our numerical
simulations. When statistically steady solutions appear, the possibility of
multifractality should be considered [22, 23]. A future line of approach might
be to repeat the calculation in [52], where the correspondence between the
multifractal model of turbulence and the NSEs was investigated.

Acknowledgments : RP and JDG would like to thank the Isaac Newton Institute
for Mathematical Sciences for support and hospitality during the programme
Mathematical Aspects of Fluid Turbulence : where do we stand? in 2022, when
work on this paper was undertaken. It was supported by grant number
EP/R014604/1. RP, KVK and NBP thank SERB (India), the National Supercomputing
Mission (India), and SERC (IISc), for support and computational resources, and
A. Gupta for discussions.

##  Appendix A Differential Inequalities

The most widely used class of differential inequalities that generalize the
Sobolev inequalities are called Gagliardo-Nirenberg inequalities [53]. In
their most general form in integer $d$ dimensions ($d=1,2,3$) they can be
expressed as

$\|\nabla^{j}u\|_{p}\leq c\,\|\nabla^{n}u\|_{r}^{a}\|u\|_{q}^{1-a}\,,$ (A.1)

where $0\leq j<n$ and $1<p,\,r,\,q\leq\infty$. The exponent $a$ can be
calculated by dimensional analysis and thus must satisfy

$\frac{1}{p}=\frac{j}{d}+a\left(\frac{1}{r}-\frac{n}{d}\right)+\frac{1-a}{q}\,,$
(A.2)

where $j/n\leq a<1$. (A.2) holds on the whole space $\mathbb{R}^{d}$. With
periodic boundary conditions, there are $L^{2}$ additional terms there are
lower-order corrections to our estimates, which will be ignored. In the
following :

  1. 1.

The standard Leray-Hopf weak-solution machinery has been used, which has been
derived for the NSEs in [31].

  2. 2.

We will use the convention that the constants designated as $c$, $c_{m}$ or
$c_{n,m}$ are generic in the sense that they may differ from line to line.

##  Appendix B Proofs in the $d=2$ case

###  B.1 Estimates for $\left<P_{2,1}\right>_{T}$

Our first requirement is to bound $P_{2,1}$, with $P_{n,1}$ defined in (3.12).
Clearly $\alpha_{n,1,2}=1/n$, so

$\displaystyle\left<P_{2,1}\right>_{T}=\left<H_{2}^{1/2}\right>_{T}$
$\displaystyle=$
$\displaystyle\left<\left(\frac{H_{2}}{H_{1}}\right)^{1/2}H_{1}^{1/2}\right>_{T}$
(B.1) $\displaystyle\leq$
$\displaystyle\left<\frac{H_{2}}{H_{1}}\right>_{T}^{1/2}\left<P_{1,1}\right>_{T}^{1/2}\,.$

To achieve a bound on $\left<H_{2}/H_{1}\right>_{T}$ , we take the curl of the
ITT equation :

$\left(\partial_{t}+\mbox{\boldmath$u$}\cdot\nabla\right)\mbox{\boldmath$\omega$}=\alpha_{0}\mbox{\boldmath$\omega$}+{\rm
Re}_{\nu}^{-1}\Delta\mbox{\boldmath$\omega$}-{\rm
Re}_{\beta}\,\mbox{curl}(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2})\,.$
(B.2)

The key point is that, while the vortex-stretching term
$\mbox{\boldmath$\omega$}\cdot\nabla\mbox{\boldmath$u$}$ is missing because of
the orthogonality of $\omega$ with the gradient operator, there is an
additional negative ${\rm
Re}_{\beta}\,\mbox{curl}(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2})$-term.
To deal with this we write

$\displaystyle{\scriptstyle\frac{1}{2}}\dot{H}_{1}$ $\displaystyle=$
$\displaystyle\alpha_{0}H_{1}-{\rm Re}_{\nu}^{-1}H_{2}-{\rm
Re}_{\beta}\int_{V}\mbox{\boldmath$\omega$}\cdot\mbox{curl}(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2})dV$
(B.3) $\displaystyle\leq$
$\displaystyle\alpha_{0}H_{1}-{\scriptstyle\frac{1}{2}}{\rm
Re}_{\nu}^{-1}H_{2}+{\scriptstyle\frac{1}{2}}{\rm Re}_{\beta}^{2}{\rm
Re}_{\nu}\int_{V}|\mbox{\boldmath$u$}|^{6}dV\,,$

where we have integrated by parts and have then used a Hölder inequality. Now
divide by $H_{1}$ to obtain

${\scriptstyle\frac{1}{2}}\left<\frac{H_{2}}{H_{1}}\right>_{T}\leq{\rm
Re}_{\nu}\alpha_{0}+{\scriptstyle\frac{1}{2}}{\rm Re}_{\beta}^{2}{\rm
Re}_{\nu}^{2}\left<\frac{\int_{V}|\mbox{\boldmath$u$}|^{6}dV}{H_{1}}\right>_{T}\,.$
(B.4)

As in (6.13), we use the Gagliardo-Nirenberg inequality

$\|\mbox{\boldmath$u$}\|_{6}\leq
c\,\|\nabla\mbox{\boldmath$u$}\|_{2}^{a}\|\mbox{\boldmath$u$}\|_{4}^{1-a}\,,\qquad
a={\scriptstyle\frac{1}{3}}\,,$ (B.5)

to find

$\left<\frac{\|\mbox{\boldmath$u$}\|_{6}^{6}}{H_{1}}\right>_{T}\leq
c\,\left<\|\mbox{\boldmath$u$}\|_{4}^{4}\right>_{T}\,.$ (B.6)

Inserting this into (B.4) and using (B.5) and the definition of
$\mathcal{A}_{0}$ in (4.9), we find that

$\displaystyle{\scriptstyle\frac{1}{2}}\left<\frac{H_{2}}{H_{1}}\right>_{T}$
$\displaystyle\leq$ $\displaystyle\alpha_{0}{\rm Re}_{\nu}+c\,{\rm
Re}_{\beta}^{2}{\rm
Re}_{\nu}^{2}\left<\|\mbox{\boldmath$u$}\|_{4}^{4}\right>_{T}$ (B.7)
$\displaystyle\leq$ $\displaystyle c\,\alpha_{0}{\rm
Re}_{\nu}\left(1+\alpha_{0}{\rm Re}_{\nu}\right)\,.$

Thus, to leading order in ${\rm Re}_{\nu}$, (B.1) becomes

$\left<P_{2,1}\right>_{T}\lessapprox
c\,\mathcal{A}_{0}^{1/2}\left(\alpha_{0}{\rm Re}_{\nu}\right)^{3/2}\,,$ (B.8)

as advertised in (6.2).

For the second option in which $U_{0}$ is chosen to be $U_{0}=\nu/L$, thus
making ${\rm Re}_{\nu}=1$, a re-working of the bounds above gives

$\left<P_{2,1}\right>_{T}\leq
c\,\alpha_{0}^{3/2}\left(1+\alpha_{0}\right)^{1/2}{\rm Re}_{\beta}^{-1/2}\,.$
(B.9)

We leave the remaining members of this second class of estimates to be
calculated by the reader.

###  B.2 An estimate for
$P_{1,m}=\left<\|\nabla\mbox{\boldmath$u$}\|_{2m}^{\frac{2m}{2m-1}}\right>_{T}$

A Gagliardo-Nirenberg inequality shows that, for some function $A$, for $m\geq
1$

$\|A\|_{2m}\leq c_{m}\|\nabla
A\|_{2}^{\frac{m-1}{m}}\|A\|_{2}^{\frac{1}{m}}\,.$ (B.10)

We choose $A=\nabla\mbox{\boldmath$u$}$ and, noting from (3.4) that
$\alpha_{1,m,2}=m/(2m-1)$, we write

$\displaystyle\left<\|\nabla\mbox{\boldmath$u$}\|_{2m}^{\frac{2m}{2m-1}}\right>_{T}$
$\displaystyle\leq$ $\displaystyle
c_{m}\left<\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}^{\frac{2(m-1)}{2m-1}}\|\nabla\mbox{\boldmath$u$}\|_{2}^{\frac{2}{2m-1}}\right>_{T}$
(B.11) $\displaystyle\leq$ $\displaystyle
c_{m}\left<\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}\right>_{T}^{\frac{2(m-1)}{2m-1}}\left<\|\nabla\mbox{\boldmath$u$}\|_{2}^{2}\right>_{T}^{\frac{1}{2m-1}}$
$\displaystyle=$ $\displaystyle
c_{m}\left<P_{2,1}\right>_{T}^{\frac{2(m-1)}{2m-1}}\left<H_{1}\right>_{T}^{\frac{1}{2m-1}}$
$\displaystyle\leq$ $\displaystyle
c_{m}\mathcal{A}_{0}^{\frac{m}{2m-1}}\left(\alpha_{0}{\rm
Re}_{\nu}\right)^{\frac{3m-2}{2m-1}}\,,$

as advertised in (6.5). In the limit $m\to\infty$ ,

$\left<\|\nabla\mbox{\boldmath$u$}\|_{\infty}\right>_{T}\leq
c\,\mathcal{A}_{0}^{1/2}\left(\alpha_{0}{\rm Re}_{\nu}\right)^{3/2}\,.$ (B.12)

###  B.3 Estimates for
$\left<P_{0,m}\right>_{T}=\left<\|\mbox{\boldmath$u$}\|_{2m}^{\frac{2m}{m-1}}\right>_{T}$
and $\left<\|\mbox{\boldmath$u$}\|_{\infty}^{2}\right>_{T}$

We now turn to estimating $u$ in $L^{2m}(V)$ for $m>2$. A Gagliardo-Nirenberg
inequality shows that

$\|\mbox{\boldmath$u$}\|_{2m}\leq
c_{m}\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}^{a}\|\mbox{\boldmath$u$}\|_{4}^{1-a}\,,$
(B.13)

where $a=\frac{m-2}{3m}$. When $n=0$ and $d=2$ we have
$(4-d)\alpha_{0,m,2}=\frac{2m}{m-1}$. Thus,

$\displaystyle\left<\|\mbox{\boldmath$u$}\|_{2m}^{\frac{2m}{m-1}}\right>_{T}$
$\displaystyle\leq$ $\displaystyle
c_{m}\left<\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}^{\frac{2(m-2)}{3(m-1)}}\|\mbox{\boldmath$u$}\|_{4}^{\frac{4(m+1)}{3(m-1)}}\right>_{T}$
(B.14) $\displaystyle\leq$
$\displaystyle\left<P_{2,1}\right>_{T}^{\frac{2(m-2)}{3(m-1)}}\left<\|\mbox{\boldmath$u$}\|_{4}^{4}\right>_{T}^{\frac{(m+1)}{3(m-1)}}\,,$

in which case, for $m>2$, using (B.8) and (4.9),

$\left<\|\mbox{\boldmath$u$}\|_{2m}^{\frac{2m}{m-1}}\right>_{T}\leq
c\,\mathcal{A}_{0}^{\frac{m}{m-1}}\left(\alpha_{0}{\rm
Re}_{\nu}\right)^{\frac{m-2}{m-1}}\,,$ (B.15)

as advertised in (6.8). In the limit $m\to\infty$ , this reduces to

$\left<\|\mbox{\boldmath$u$}\|_{\infty}^{2}\right>_{T}\leq
c\,\alpha_{0}\mathcal{A}_{0}{\rm Re}_{\nu}\,,$ (B.16)

as advertised in (6.9).

###  B.4 Estimates for $\left<P_{n,1}\right>_{T}$ for $n>2$

Using the methods in [29], a full ‘ladder’ for $H_{n}$ takes the form

$\displaystyle{\scriptstyle\frac{1}{2}}\dot{H}_{n}$ $\displaystyle\leq$
$\displaystyle\alpha_{0}H_{n}-{\rm
Re}_{\nu}^{-1}H_{n+1}+c_{n,1}H_{n+1}^{1/2}H_{n}^{1/2}\|\mbox{\boldmath$u$}\|_{\infty}$
(B.17) $\displaystyle+$ $\displaystyle c_{n,2}{\rm
Re}_{\beta}H_{n}\|\mbox{\boldmath$u$}\|_{\infty}^{2}\,;$

therefore, after a Hölder inequality and re-arrangement, we have

${\scriptstyle\frac{1}{2}}\dot{H}_{n}\leq\alpha_{0}H_{n}-{\scriptstyle\frac{1}{2}}{\rm
Re}_{\nu}^{-1}H_{n+1}+c_{n}\left({\rm Re}_{\beta}+{\rm
Re}_{\nu}\right)H_{n}\|\mbox{\boldmath$u$}\|_{\infty}^{2}\,.$ (B.18)

Thus,

$\left<\frac{H_{n+1}}{H_{n}}\right>_{T}\leq 2\alpha_{0}{\rm
Re}_{\nu}+c_{n}{\rm Re}_{\nu}\left({\rm Re}_{\beta}+{\rm
Re}_{\nu}\right)\left<\|\mbox{\boldmath$u$}\|_{\infty}^{2}\right>_{T}\,.$
(B.19)

From this and (B.16) we deduce that

$\displaystyle\left<\frac{H_{n+1}}{H_{n}}\right>_{T}$
$\displaystyle\lessapprox$ $\displaystyle c_{n}\alpha_{0}\mathcal{A}_{0}{\rm
Re}_{\nu}^{3}\,.$ (B.20)

Moreover, for $n\geq 2$

$\displaystyle\left<P_{n+1,1}\right>_{T}$ $\displaystyle=$
$\displaystyle\left<\left(\frac{H_{n+1}}{H_{n}}\right)^{\frac{1}{n+1}}H_{n}^{\frac{1}{n+1}}\right>_{T}$
(B.21) $\displaystyle\leq$
$\displaystyle\left<\frac{H_{n+1}}{H_{n}}\right>_{T}^{\frac{1}{n+1}}\left<P_{n,1}\right>_{T}^{\frac{n}{n+1}}\,.$

Given that $\left<P_{2,1}\right>_{T}$ is bounded above as in (B.8), then
together with (B.20), we can now generate upper bounds on every
$\left<P_{n,1}\right>_{T}$. For $n\geq 2$, these are

$\displaystyle\left<P_{n,1}\right>_{T}$ $\displaystyle\leq$ $\displaystyle
c_{n,1}\alpha_{0}^{\frac{n+1}{n}}\mathcal{A}_{0}^{\frac{n-1}{n}}{\rm
Re}_{\nu}^{\frac{3(n-1)}{n}}\,,$ (B.22)

which can be transformed into the form advertised in (6.3).

###  B.5 Estimates for $\left<P_{n,m}\right>_{T}$

We can write down an inequality for $B=\nabla^{2}\mbox{\boldmath$u$}$ such
that

$\|\nabla^{n-2}B\|_{2m}\leq c\,\|\nabla^{N-2}B\|_{2}^{a}\|B\|_{2}^{1-a}\,,$
(B.23)

for some $N>n+1-1/m$ , with

$a=\frac{m(n-1)-1}{m(n+1)-1}\,.$ (B.24)

Thus, we can write

$\left<\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}^{2\alpha_{n,m}}\right>_{T}\leq
c\,\left<\|\nabla^{N}\mbox{\boldmath$u$}\|_{2}^{2a\alpha_{n,m}}\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}^{2(1-a)\alpha_{n,m}}\right>_{T}\,.$
(B.25)

Re-arranging and then using Hölder’s inequality, we have

$\displaystyle\left<\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}^{2\alpha_{n,m}}\right>_{T}$
$\displaystyle\leq$ $\displaystyle
c\,\left<\left(\|\nabla^{N}\mbox{\boldmath$u$}\|_{2}^{2/N}\right)^{aN\alpha_{n,m}}\left(\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}\right)^{2(1-a)\alpha_{n,m}}\right>_{T}$
(B.26) $\displaystyle\leq$ $\displaystyle
c_{N,n,m}\left<\|\nabla^{N}\mbox{\boldmath$u$}\|_{2}^{2/N}\right>_{T}^{aN\alpha_{n,m}}\left<\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}^{\frac{2(1-a)\alpha_{n,m}}{1-aN\alpha_{n,m}}}\right>_{T}^{1-aN\alpha_{n,m}}$
(B.27) $\displaystyle=$ $\displaystyle
c_{N,n,m}\left<P_{N,1}\right>_{T}^{aN\alpha_{n,m}}\left<P_{2,1}^{\frac{2(1-a)\alpha_{n,m}}{1-aN\alpha_{n,m}}}\right>_{T}^{1-aN\alpha_{n,m}}\,.$
(B.28)

It can be checked through the definition of $a$ in (B.24) that the exponent of
$P_{2,1}$ inside the time-average is unity. Estimates for
$\left<P_{N,1}\right>_{T}$ and $\left<P_{2,1}\right>_{T}$ are available from
(6.3) and (6.2) : one can then choose the lowest value of $N$ , constrained by
$N>n+1-1/m$. After some algebra, this leads to the result

$\left<P_{n,m}\right>_{T}\leq
c_{n,m}\alpha_{0}^{\frac{2m}{m(n+1)-1}}\left(\alpha_{0}\mathcal{A}_{0}{\rm
Re}_{\nu}^{3}\right)^{\frac{mn-1}{m(n+1)-1}}\,,$ (B.29)

as advertised in (6.10).

##  Appendix C Proofs in the $d=3$ case

Step 1 : $Q_{n,m}$ is defined in (3.15). In simplied form $Q_{n,1}$ can be
written as

$Q_{n,1}=H_{n}^{\frac{1}{2n-1}}\,.$ (C.1)

Moreover, because $\left<Q_{1,1}\right>_{T}=\left<H_{1}\right>_{T}$ , we have
an estimate for this in (4.8). We begin this section by estimating
$\left<Q_{2,1}\right>_{T}$ from the vorticity equation (B.2), with the $3d$
vortex stretching term restored :

$\left(\partial_{t}+\mbox{\boldmath$u$}\cdot\nabla\right)\mbox{\boldmath$\omega$}=\alpha_{0}\mbox{\boldmath$\omega$}+{\rm
Re}_{\nu}^{-1}\Delta\mbox{\boldmath$\omega$}+\mbox{\boldmath$\omega$}\cdot\nabla\mbox{\boldmath$u$}-{\rm
Re}_{\beta}\,\mbox{curl}(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2})\,.$
(C.2)

The equivalent of (B.3) is

$\displaystyle{\scriptstyle\frac{1}{2}}\dot{H}_{1}$ $\displaystyle\leq$
$\displaystyle\alpha_{0}H_{1}-{\rm
Re}_{\nu}^{-1}H_{2}+c\,H_{2}^{1/2}H_{1}^{1/2}\|\mbox{\boldmath$u$}\|_{\infty}-{\rm
Re}_{\beta}\int_{V}\mbox{\boldmath$\omega$}\cdot\mbox{curl}(\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2})dV$
(C.3) $\displaystyle\leq$
$\displaystyle\alpha_{0}H_{1}-{\scriptstyle\frac{3}{4}}{\rm
Re}_{\nu}^{-1}H_{2}+c\,H_{2}^{1/2}H_{1}^{1/2}\|\mbox{\boldmath$u$}\|_{\infty}+{\rm
Re}_{\beta}^{2}{\rm Re}_{\nu}\int_{V}|\mbox{\boldmath$u$}|^{6}dV\,,$

where we have integrated by parts and have then used a Hölder inequality. The
$3$-dimensional Agmon inequality for $\|\mbox{\boldmath$u$}\|_{\infty}$ is

$\|\mbox{\boldmath$u$}\|_{\infty}^{2}\leq c_{n}H_{n}^{a}H_{1}^{1-a}\qquad
n\geq 2\,,$ (C.4)

with $a=\frac{1}{2(n-1)}$. Thus, for $n=2$ ,

$\|\mbox{\boldmath$u$}\|_{\infty}\leq c\,H_{2}^{1/4}H_{1}^{1/4}$ (C.5)

and so

$\displaystyle H_{2}^{1/2}H_{1}^{1/2}\|\mbox{\boldmath$u$}\|_{\infty}$
$\displaystyle\leq$ $\displaystyle c\,H_{2}^{3/4}H_{1}^{3/4}$ (C.6)
$\displaystyle\leq$ $\displaystyle{\scriptstyle\frac{1}{4}}{\rm
Re}_{\nu}^{-1}H_{2}+{{\scriptstyle\frac{3}{4}}}c\,{\rm
Re}_{\nu}^{3}H_{1}^{3}\,.$ (C.7)

Moreover, Sobolev’s inequality for $d=3$ shows that

$\|\mbox{\boldmath$u$}\|_{6}\leq c\,\|\nabla\mbox{\boldmath$u$}\|_{2}\,.$
(C.8)

Therefore, in total, (C.3) becomes

${\scriptstyle\frac{1}{2}}\dot{H}_{1}\leq\alpha_{0}H_{1}-{\scriptstyle\frac{1}{2}}{\rm
Re}_{\nu}^{-1}H_{2}+c\,{\rm Re}_{\nu}\left({\rm Re}_{\beta}^{2}+{\rm
Re}_{\nu}^{2}\right)H_{1}^{3}\,.$ (C.9)

Thus, the ultimate contribution to (C.10) from the
$\mbox{\boldmath$u$}|\mbox{\boldmath$u$}|^{2}$-term is proportional to that
from the vortex-stretching term, in the sense that they are both proportional
to $H_{1}^{3}$. Dividing by $H_{1}^{2}$ gives

${\frac{1}{2}}{\rm
Re}_{\nu}^{-1}\left<\frac{H_{2}}{H_{1}^{2}}\right>_{T}\leq\alpha_{0}\left<H_{1}^{-1}\right>_{T}+c\,{\rm
Re}_{\nu}\left({\rm Re}_{\beta}^{2}+{\rm
Re}_{\nu}^{2}\right)\left<Q_{1,1}\right>_{T}\,.$ (C.10)

Ignoring the first term on the right hand side with the negative exponent, we
can write

$\left<\frac{H_{2}}{H_{1}^{2}}\right>_{T}\leq c\,\alpha_{0}{\rm
Re}_{\nu}^{2}\left({\rm Re}_{\beta}^{2}+{\rm Re}_{\nu}^{2}\right)\,.$ (C.11)

Thus, we finally have

$\displaystyle\left<Q_{2,1}\right>_{T}$ $\displaystyle=$
$\displaystyle\left<\left(\frac{H_{2}}{H_{1}^{2}}\right)^{1/3}H_{1}^{2/3}\right>_{T}$
(C.12) $\displaystyle\leq$
$\displaystyle\left<\frac{H_{2}}{H_{1}^{2}}\right>_{T}^{1/3}\left<H_{1}\right>_{T}^{2/3}$
(C.13) $\displaystyle\leq$ $\displaystyle c\,\left(\alpha_{0}{\rm
Re}_{\nu}^{2}\left({\rm Re}_{\beta}^{2}+{\rm
Re}_{\nu}^{2}\right)\right)^{1/3}(\alpha_{0}{\rm Re}_{\nu})^{2/3}$ (C.14)
$\displaystyle=$ $\displaystyle c\,\alpha_{0}{\rm Re}_{\nu}^{4/3}\left({\rm
Re}_{\beta}^{2}+{\rm Re}_{\nu}^{2}\right)^{1/3}\,.$ (C.15)

If ${\rm Re}_{\nu}$ is dominant, the bound scales like

$\left<Q_{2,1}\right>_{T}\lessapprox\alpha_{0}{\rm Re}_{\nu}^{2}+O\left({\rm
Re}_{\nu}^{4/3}\right)\,.$ (C.16)

Step 2 : Let us repeat (B.18) by writing

${\scriptstyle\frac{1}{2}}\dot{H}_{n}\leq\alpha_{0}H_{n}-{\scriptstyle\frac{1}{2}}{\rm
Re}_{\nu}^{-1}H_{n+1}+c_{n}\left({\rm Re}_{\beta}+{\rm
Re}_{\nu}\right)H_{n}\|\mbox{\boldmath$u$}\|_{\infty}^{2}\,.$ (C.17)

After re-arrangement and the use of Agmon’s inequality, (C.17) becomes

$\displaystyle{\scriptstyle\frac{1}{2}}\dot{H}_{n}\leq\alpha_{0}H_{n}-{\scriptstyle\frac{1}{2}}{\rm
Re}_{\nu}^{-1}H_{n+1}+c_{n}\left({\rm Re}_{\beta}+{\rm
Re}_{\nu}\right)H_{n}^{1+a}H_{1}^{1-a}\,.$ (C.18)

Dividing by $H_{n}^{2n/(2n-1)}$ and time averaging gives

$\displaystyle{\scriptstyle\frac{1}{2}}{\rm
Re}_{\nu}^{-1}\left<\frac{H_{n+1}}{H_{n}^{\frac{2n}{2n-1}}}\right>_{T}$
$\displaystyle\leq$
$\displaystyle\alpha_{0}\left<H_{n}^{1-\frac{2n}{2n-1}}\right>_{T}+c_{n}\left({\rm
Re}_{\beta}+{\rm
Re}_{\nu}\right)\left<H_{n}^{\frac{2n-1}{2(n-1)}-\frac{2n}{2n-1}}H_{1}^{\frac{2n-3}{2(n-1)}}\right>_{T}$
(C.19) $\displaystyle\leq$
$\displaystyle\alpha_{0}\left<H_{n}^{-\frac{1}{2n-1}}\right>_{T}+c_{n}\left({\rm
Re}_{\beta}+{\rm
Re}_{\nu}\right)\left<Q_{n}\right>_{T}^{\frac{1}{2(n-1)}}\left<Q_{1}\right>^{\frac{2n-3}{2(n-1)}}\,.$

The next step is to ignore the first term555The term with the negative
exponent of $H_{n}$ on the RHS of (C.19) is only out of control if $H_{n}$
temporarily becomes very small. In principle, this could be dealt with by
adding a constant term to $H_{n}$ to provide the platform of a lower bound. We
omit the details.. Given that ${\rm Re}_{\nu}$ is the dominant term, we write
(C.19) in the simplified form

$\left<\frac{H_{n+1}}{H_{n}^{\frac{2n}{2n-1}}}\right>_{T}\leq c_{n}{\rm
Re}_{\nu}^{2}\left<Q_{n,1}\right>_{T}^{\frac{1}{2(n-1)}}\left<Q_{1,1}\right>_{T}^{\frac{2n-3}{2(n-1)}}\,.$
(C.20)

We then study

$\displaystyle\left<Q_{n+1,1}\right>_{T}=\left<H_{n+1}^{\frac{1}{2n+1}}\right>_{T}$
$\displaystyle=$
$\displaystyle\left<\left(\frac{H_{n+1}}{H_{n}^{\frac{2n}{2n-1}}}\right)^{\frac{1}{2n+1}}H_{n}^{\frac{2n}{(2n+1)(2n-1)}}\right>_{T}$
(C.21) $\displaystyle\leq$
$\displaystyle\left<\frac{H_{n+1}}{H_{n}^{\frac{2n}{2n-1}}}\right>_{T}^{\frac{1}{2n+1}}\left<Q_{n,1}\right>_{T}^{\frac{2n}{2n+1}}\,,$

in which case

$\displaystyle\left<Q_{n+1,1}\right>_{T}$ $\displaystyle\leq$ $\displaystyle
c_{n}{\rm
Re}_{\nu}^{\frac{2}{2n+1}}\left<Q_{n,1}\right>_{T}^{\frac{2n}{2n+1}+\frac{1}{2(2n+1)(n-1)}}\left<Q_{1,1}\right>_{T}^{\frac{2n-3}{2(n-1)(2n+1)}}$
(C.22) $\displaystyle=$ $\displaystyle c_{n}{\rm
Re}_{\nu}^{\frac{2}{2n+1}}\left<Q_{n,1}\right>_{T}^{\frac{(2n-1)^{2}}{2(2n+1)(n-1)}}\left<Q_{1,1}\right>_{T}^{\frac{2n-3}{2(n-1)(2n+1)}}\,.$

Given that we have estimates for both $\left<Q_{1,1}\right>_{T}$ and
$\left<Q_{2,1}\right>_{T}$ , we can generate estimates for all
$\left<Q_{n,1}\right>_{T}$ for $n\geq 3$. Thus, we deduce that

$\left<Q_{n,1}\right>_{T}<\infty\qquad n\geq 3\,.$ (C.23)

The bound is messy so we simply register that the right hand side is finite.
Finally, the method used in §B.5 can be used to show that
$\left<Q_{n,m}\right>_{T}<\infty$ for $m\geq 1$.

Step 3 : Now let us consider

$\|\mbox{\boldmath$u$}\|_{2m}\leq
c\,\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}^{a}\|\mbox{\boldmath$u$}\|_{4}^{1-a}\,,$
(C.24)

where $a=3(m-2)/5m$ with $m>2$. Because
$\alpha_{2,1,3}={\scriptstyle\frac{2}{3}}$ we can write

$\displaystyle\left<\|\mbox{\boldmath$u$}\|_{2m}^{\alpha_{0,m}}\right>_{T}$
$\displaystyle\leq$ $\displaystyle
c\,\left<\|\nabla^{2}\mbox{\boldmath$u$}\|_{2}^{a\alpha_{0,m}}\|\mbox{\boldmath$u$}\|_{4}^{(1-a)\alpha_{0,m}}\right>_{T}$
(C.25) $\displaystyle=$ $\displaystyle
c\,\left<Q_{2,1}^{3a\alpha_{0,m}/2}\left(\|\mbox{\boldmath$u$}\|_{4}^{4}\right)^{{\scriptstyle\frac{1}{4}}(1-a)\alpha_{0,m}}\right>_{T}\,,$

where $\alpha_{0,m}=\frac{2m}{2m-3}$. Then

$\left<\|\mbox{\boldmath$u$}\|_{2m}^{\alpha_{0,m}}\right>_{T}\leq
c\,\left<Q_{2,1}\right>_{T}^{3a\alpha_{0,m}/2}\left<\left(\|\mbox{\boldmath$u$}\|_{4}^{4}\right)^{\frac{{\scriptstyle\frac{1}{4}}(1-a)\alpha_{0,m}}{1-3a\alpha_{0,m}/2}}\right>_{T}^{1-3a\alpha_{0,m}/2}\,.$
(C.26)

Given $a$ and $\alpha_{0,m}$, it can easily be checked that the exponent of
$\|\mbox{\boldmath$u$}\|_{4}^{4}$ inside the average is unity. Thus, because
$\left<\|\mbox{\boldmath$u$}\|_{4}^{4}\right>_{T}\leq c\,\mathcal{A}_{0}^{2}$
and with the help of (C.16), (C.26) becomes

$\displaystyle\left<\|\mbox{\boldmath$u$}\|_{2m}^{\alpha_{0,m}}\right>_{T}$
$\displaystyle\leq$ $\displaystyle c\,(\alpha_{0}{\rm
Re}_{\nu}^{2})^{3a\alpha_{0,m}/2}\mathcal{A}_{0}^{2-3a\alpha_{0,m}}$ (C.27)
$\displaystyle=$
$\displaystyle\alpha_{0}^{3a\alpha_{0,m}/2}\mathcal{A}_{0}^{2-3a\alpha_{0,m}}{\rm
Re}_{\nu}^{3a\alpha_{0,m}}\,.$

In fact $3a\alpha_{0,m}/2=\frac{9(m-2)}{5(2m-3)}$ and so
$1-3a\alpha_{0,m}/2=\frac{m+3}{5(2m-3)}$, whence

$\left<\|\mbox{\boldmath$u$}\|_{2m}^{\alpha_{0,m}}\right>_{T}\leq
c\,\mathcal{A}_{0}^{\frac{2(m+3)}{5(2m-3)}}\left(\alpha_{0}{\rm
Re}_{\nu}^{2}\right)^{\frac{9(m-2)}{5(2m-3)}}\qquad m>2\,,$ (C.28)

as advertised in (7.4). In the limit $m\to\infty$ , we find that

$\displaystyle\left<\|\mbox{\boldmath$u$}\|_{\infty}\right>_{T}$
$\displaystyle\leq$ $\displaystyle
c\,\mathcal{A}_{0}^{1/5}\left(\alpha_{0}{\rm Re}_{\nu}^{2}\right)^{9/10}\,,$
(C.29)

as advertised in (7.5).

## References

  * [1]
  * [2] C. W. Reynolds, Flocks, Herds, and Schools : A Distributed Behavioral Model, Computer Graphics, 21, No. 4 (1987) pp. 25–34.
  * [3] T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel Type of Phase Transition in a System of Self-Driven Particles, Phys. Rev. Lett. 75 (1995) 1226–1229.
  * [4] A. Czirok, H. E. Stanley, and T. Vicsek, Spontaneously ordered motion of self-propelled particles, J. Phys. A 30 (1997) 1375–1385.
  * [5] J. Toner and Y. Tu, Long-Range Order in a Two-Dimensional Dynamical XY-Model : How Birds Fly Together, Phys. Rev. Lett. 75 (1995) 4326–4329.
  * [6] J. Toner, Y. Tu, and S. Ramaswamy, Hydrodynamics and phases of flocks, Ann. Phys. 318 (2005) 170–244.
  * [7] R. Alert, J.-F. Joanny, and J. Casademunt, Universal scaling of active nematic turbulence, Nature Physics 16 (2020) 682–688. doi.org/10.1038/s41567-020-0854-4 .
  * [8] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Meso-scale turbulence in living fluids, Proc. Nat. Acad. Sci. 109 (2012) 14308–13.
  * [9] J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär, and R. E. Goldstein, Fluid dynamics of bacterial turbulence, Phys. Rev. Lett. 110 (2013) 228102–05.
  * [10] J. Slomka and J. Dunkel, Generalized Navier-Stokes equations for active suspensions, Euro Phys. J. Special Topics, 224 (2015) 1349–1358.
  * [11] M. Linkmann, G. Boffetta, M. C. Marchetti, and B. Eckhardt, Phase transition to large scale coherent structures in two-dimensional active matter turbulence, Phys. Rev. Lett. 122 (2019) 214503.
  * [12] M. Linkmann, M. C. Marchetti, G. Boffetta, and B. Eckhardt, Condensate formation and multiscale dynamics in two-dimensional active suspensions, Phys. Rev. E 101 (2020) 022609–16.
  * [13] J. Slomka and J. Dunkel, Spontaneous mirror-symmetry breaking induces inverse energy cascade in 3d active fluids, Proc. Nat. Acad. Sci. 114 (2017) 2119–2121.
  * [14] V. Bratanov, F. Jenko, and E. Frey, New class of turbulence in active fluids, Proc. Nat. Acad. Sci. 112 (2015) 15048–15053.
  * [15] R. Alert, J. Casademunt, and J.-F. Joanny, Active turbulence, Ann. Rev. Cond. Matter Phys. 13 (2021) 143-170.
  * [16] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. Liverpool, and J. Prost, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85 (2013) 1143–1189.
  * [17] L. Giomi, Geometry and topology of turbulence in active nematics, Phys. Rev. X 5 (2015) 031003–9.
  * [18] C. P. Sanjay and A. Joy, Friction scaling laws for transport in active turbulence, Phys. Rev. Fluids, 5 (2020) 024302.
  * [19] S. Mukherjee, R. K. Singh, M. James, S. S. Ray, Anomalous diffusion and Lévy walks distinguish active from inertial turbulence, Phys. Rev. Lett. 127 (2021) 118001–4.
  * [20] K. V. Kiran, A. Gupta, A. K. Verma, R. Pandit, Irreversiblity in Bacterial Turbulence: Insights from the Mean-Bacterial-Velocity Model (2022) arXiv:2201.12722.
  * [21] J. Swift, P. C. Hohenberg, Hydrodynamic fluctuations at the convective instability, Phys. Rev. A 15 (1977) 319–328 ; M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys., 65, No. 3, (1993).
  * [22] N. Rana, P. Perlekar, Coarsening in the two-dimensional incompressible Toner-Tu equation : Signatures of turbulence, Phys. Rev. E 102 (2020) 032617-9.
  * [23] N. Rana, P. Perlekar, Phase ordering, topological defects, and turbulence in the three-dimensional incompressible Toner-Tu equation, Phys. Rev. E 105 (2022) L032603-6.
  * [24] L. Chen, C.F. Lee, J. Toner, Mapping two-dimensional polar active ﬂuids to two-dimensional soap and one-dimensional sandblasting, Nature Comm. 7 (2016) 12215-10.
  * [25] L. Chen, C. F. Lee, J. Toner, Incompressible polar active ﬂuids in the moving phase in dimensions $d>2$, New J. Phys. 20 (2018) 113035.
  * [26] M. J. Bowick, N. Fakhri, M. C. Marchetti, S. Ramaswamy, Symmetry, Thermodynamics, and Topology in Active Matter, Phys. Rev. X 12 (2022) 010501-27.
  * [27] J. Leray, Sur le mouvement d’un liquide visqueux emplissant l’espace, Acta Math., 63 (1934) 193–248.
  * [28] C. Foias, O. Manley, R. Rosa, R. Temam, Navier-Stokes Equations and Turbulence, Cambridge University Press (Cambridge 2001).
  * [29] C. Doering, J. D. Gibbon, Applied Analysis of the Navier-Stokes equations, Cambridge University Press (Cambridge 1995).
  * [30] J. C. Robinson, J. L. Rodrigo and W. Sadowski, The Three-dimensional Navier-Stokes Equations : Classical Theory ; Cambridge Studies in Advanced Mathematics. Cambridge University Press, (Cambridge 2016).
  * [31] C. Foias, C. Guillopé, and R. Temam, New a priori estimates for Navier-Stokes equations in dimension 3, Comm. PDEs 6 (1981) 329–359.
  * [32] J. D. Gibbon, Weak and strong solutions of the 3D Navier-Stokes equations and their relation to a chessboard of convergent inverse length scales, J. Nonlin. Sci. 29 (2019) 215–228.
  * [33] J. D. Gibbon, Turbulent cascades and thin sets in 3D Navier-Stokes turbulence, Europhys. Lett. 131 (2020) 64001–05.
  * [34] F. Zanger, H. Löwen, J. Saal, MNP2015 2015: Math. Nonlin. Phen. – Anal. Comp. (2017) 285–-303.
  * [35] S. Antontsev, H. de Oliveira, The Navier–Stokes problem modified by an absorption term, App. Anal. 89 (2010) 1805–1825.
  * [36] E. S. Titi, S. Trabelsi, Global well-posedness of a three-dimensional Brinkman-Forchheimer-Bénard convection model in porous media, (2022) arXiv:2204.03531v1.
  * [37] X. Cai, Q. Jiu, Weak and strong solutions for the incompressible Navier–Stokes equations with damping, J. Math. Anal. App. 343 (2008) 799–809.
  * [38] A. O. Celebi, V. K. Kalantarov, D. Ugurlu, On continuous dependence on coefficients of the Brinkman–Forchheimer equations, App. Math. Lett. 19 (2006) 801–807.
  * [39] Y. Liu, C. Lin, Structural stability for Brinkman–Forchheimer equations., Elec. J. Diff. Equs. 2007(02), (2007) 1–8.
  * [40] M. Louaked, N. Seloula, S. Sun, S. Trabelsi, A pseudocompressibility method for the incompressible Brinkman-Forchheimer equations, Differ. Integral Equ. 28 (2015) 361–382.
  * [41] M. Louaked, N. Seloula, S. Trabelsi, Approximation of the unsteady Brinkman-Forchheimer equations by the pressure stabilization method, Numer. Meth. PDEs 33 (2017) 1949–1965.
  * [42] P. A. Markowich, E. S. Titi, S. Trabelsi, Continuous data assimilation for the three-dimensional Brinkman–Forchheimer-extended Darcy model, Nonlinearity 29 (2016) 1292–1328.
  * [43] L. E. Payne, B. Straughan, Convergence and continuous dependence for the Brinkman–Forchheimer equations, Stud. App. Math. 102 (1999) 419–439.
  * [44] B. Straughan, Stability And Wave Motion In Porous Media, 165, Springer Science & Business Media (2008).
  * [45] D. Ugurlu, On the existence of a global attractor for the Brinkman–Forchheimer equations, Nonlin. Anal : Theory, Methods & Applications 68 (2008) 1986–1992.
  * [46] B. Wang, S. Lin, Existence of global attractors for the three-dimensional Brinkman–Forchheimer equation, Math. Meth. App. Sci. 31 (2008) 1479–1495.
  * [47] Y. You, C. Zhao, S. Zhou, The existence of uniform attractors for 3d Brinkman–Forchheimer equations, Dis. Cont. Dyn. Sys. A 32 (2012) 787–3800.
  * [48] D. Donzis, R. M. Kerr, J. D. Gibbon, A. Gupta, R. Pandit, D. Vincenzi, Vorticity moments in four numerical simulations of the 3D Navier-Stokes equations, J. Fluid Mech. 732 (2013) 316–331.
  * [49] J. D. Gibbon, D. Donzis, R. M. Kerr, A. Gupta, R. Pandit, D. Vincenzi, Regimes of nonlinear depletion and regularity in the 3D Navier-Stokes equations, Nonlinearity 27 (2014) 1–19.
  * [50] J. D. Gibbon, A. Gupta, G. Krstulovic, R. Pandit, H. Politano, Y. Ponty, A. Pouquet, G. Sahoo, J. Stawarz, Depletion of Nonlinearity in Magnetohydrodynamic Turbulence: Insights from Analysis and Simulations, Phys. Rev. E 93 (2016) 043104.
  * [51] J. D. Gibbon, A. Gupta, N. Pal, R. Pandit, The role of BKM-type theorems in 3d Euler, Navier-Stokes and Cahn-Hilliard analysis, Physica D : Nonlinear Phenomena 376 (2018) 60–68.
  * [52] B. Dubrulle, J. D. Gibbon, A correspondence between the multifractal model of turbulence and the Navier–Stokes equations, Phil. Trans. R. Soc. A, 380 (2022) 20210092.
  * [53] R. A. Adams, Sobolev Spaces, Pure and Applied Mathematics, 65, Academic Press (NY, 1975).
  * [54] C. Canuto, M. Y. Hussaini, M. Yousuff, A. Quarteroni, A. Thomas Jr., Spectral methods in fluid dynamics, Springer Science & Business Media, (2012).
  * [55] S. M. Cox, P. C. Matthews, Exponential time differencing for stiff systems, J. Comput. Phys., 176 (2002) 430–455.

In this Supplemental Material we provide the following additional plots:

  1. 1.

Different norms and their time averages for the temporally frozen states in
$d=2$.

  2. 2.

In $d=3$, we plot versus $t$, $Q_{0,m}$, $Q_{1,m}$, and $D_{m}$ and their time
averages versus $m$, for a representative run. Weighted averages of the
vorticity are defined as

$D_{m}(t)=\|\mbox{\boldmath$\omega$}\|^{\alpha_{n,m,d=3}}_{2m}\,.$

##  Appendix D Supplemental results for two dimensions

AB Figure 8: (Colour online) Illustrative plots for
$U_{0}=\sqrt{\alpha/\beta}$ (panel A) and $U_{0}=\nu/L$ (panel B) for run F7
(see Table 1, main text): First and second rows : plots versus $t$ of
$P_{0,m}$ and $P_{1,m}$ ; the plots in the second row are expanded versions of
small segments of the plots in the first row. Third row : Plots versus $m$ of
$\left<P_{0,m}\right>_{T}$ and $\left<P_{1,m}\right>_{T}$. Curves for
$m=2,3,4,5,6,7,8,9$, and $10$ are drawn in red, pink; violet, green, cyan,
maroon, blue, orange, and yellow, respectively.

AB Figure 9: (Colour online) Illustrative plots for
$U_{0}=\sqrt{\alpha/\beta}$ (panel A) and $U_{0}=\nu/L$ (panel B) for $d=2$,
runs F1-F7 (see Table 1, main text) : First row : plots versus $Re_{\nu}$
(panel A) and $\alpha_{0}$ (panel B) of $\left<P_{1,1}\right>_{T}$ ( solid
black line) and $Re_{\nu}\,\alpha_{0}\,\mathcal{A}_{0}$ (dashed black line).
Second row : Plots versus $Re_{\nu}$ (panel A) and $\alpha_{0}$ (panel B) of
$\left<P_{0,m}\right>_{T}$ and $\left<P_{1,m}\right>_{T}$. Curves for
$m=2,3,4,5,6,7,8,9$, and $10$ are drawn in red, pink; violet, green, cyan,
maroon, blue, orange, and yellow, respectively.

##  Appendix E Supplemental results for three dimensions

A Figure 10: (Colour online) Illustrative plots for
$U_{0}=\sqrt{\alpha/\beta}$ (panel A) run B3 (see Table 1, main text): First
and second rows : plots versus $t$ of $Q_{0,m}$, $Q_{1,m}$ and $D_{m}$ ; the
plots in the second row are expanded versions of small segments of the plots
in the first row. Third row : Plots versus $m$ of $\left<Q_{0,m}\right>_{T}$
and $\left<Q_{1,m}\right>_{T}$ and $\left<D_{m}\right>_{T}$. Curves for
$m=2,3,4,5,6,7,8,9$, and $10$ are drawn in red, pink; violet, green, cyan,
maroon, blue, orange, and yellow, respectively. B Figure 11: (Colour online)
Illustrative plots for $U_{0}=\nu/L$ (panel B) for run B3 (see Table 1, main
text): First and second rows : plots versus $t$ of $Q_{0,m}$, $Q_{1,m}$ and
$D_{m}$ ; the plots in the second row are expanded versions of small segments
of the plots in the first row. Third row : Plots versus $m$ of
$\left<Q_{0,m}\right>_{T}$ and $\left<Q_{1,m}\right>_{T}$ and
$\left<D_{m}\right>_{T}$. Curves for $m=2,3,4,5,6,7,8,9$, and $10$ are drawn
in red, pink; violet, green, cyan, maroon, blue, orange, and yellow,
respectively.


# A novel cryogenic VUV spectrofluorometer for the characterization of
wavelength shifters

A. Leonhardt,11footnotetext: Corresponding author M. Goldbrunner  B. Hackett
and S. Schönert

###### Abstract

We present a novel cryogenic VUV spectrofluorometer designed for the
characterization of wavelength shifters (WLS) crucial for experiments based on
liquid argon (LAr) scintillation light detection. Wavelength shifters like
1,1,4,4-tetraphenyl-1,3-butadiene (TPB) or polyethylene naphthalate (PEN) are
used in these experiments to shift the VUV scintillation light to the visible
region. Precise knowledge of the optical properties of the WLS at liquid
argon’s temperature ($87\text{\,}\mathrm{K}$) and LAr scintillation wavelength
($128\text{\,}\mathrm{nm}$) is necessary to model and understand the detector
response. The cryogenic VUV spectrofluorometer was commissioned to measure the
emission spectra and relative wavelength shifting efficiency (WLSE) of samples
between $300\text{\,}\mathrm{K}87\text{\,}\mathrm{K}$ for VUV
($120\text{\,}\mathrm{nm}190\text{\,}\mathrm{nm}$) and UV
($310\text{\,}\mathrm{nm}$) excitation. New mitigation techniques for surface
effects on cold WLS were established. As part of this work, the TPB-based
wavelength shifting reflector (WLSR) featured in the neutrinoless double-beta
decay experiment LEGEND-200 was characterized. The wavelength shifting
efficiency was observed to increase by $54(4)\text{\,}\mathrm{\char 37\relax}$
from room temperature (RT) to $87\text{\,}\mathrm{K}$. PEN installed in
LEGEND-200 was also characterized and a first measurement of the relative
wavelength shifting efficiency and emission spectrum at RT and
$87\text{\,}\mathrm{K}$ is presented. Surface effects from cooling were
corrected by normalizing the measurement at VUV excitation with the respective
ones at UV excitation. The WLSE of amorphous PEN was found to be enhanced by
$49(1)\text{\,}\mathrm{\char 37\relax}$ at $87\text{\,}\mathrm{K}$ compared to
RT.

##  1 Introduction

Liquid argon (LAr) is a common scintillating medium in dark matter WIMP
detectors, such as DEAP-3600[1] and DarkSide-20k[2], and in neutrinoless
double-beta decay experiments like LEGEND-200[3] for background suppression.
Ionizing particles in LAr induce vacuum-ultraviolet (VUV) scintillation light
at $128\text{\,}\mathrm{nm}$. Most light detectors have low sensitivity in the
VUV range, necessitating the use of wavelength shifters (WLS). The widely used
1,1,4,4-tetraphenyl-1,3-butadiene (TPB) is applied as a coating to optically
inactive surfaces or to photo-detectors like photo-multiplier tubes. An
alternative is the scintillating and wavelength shifting polymer polyethylene
naphthalate (PEN), increasingly popular for its commercial availability as
pellets, which can be injection molded into structural parts or extruded into
thin films[4].

The LEGEND-200 experiment investigates neutrinoless double-beta decay of 76Ge
using $200\text{\,}\mathrm{kg}$ of enriched high-purity germanium (HPGe)
crystals, which serve as both the decay source and detectors[3]. A LAr volume,
acting as a coolant and active shield, surrounds the HPGe detectors. The light
collection system consists of light-guiding fibers surrounding the detectors,
coated in TPB and coupled to silicon photomultipliers (SiPMs). The LAr
detector’s active volume is confined by a $1.4\text{\,}\mathrm{m}$ diameter,
$3\text{\,}\mathrm{m}$ height wavelength shifting reflector (WLSR) consisting
of a Tetratex® lining coated in-situ with TPB[5]. The WLSR shifts LAr
scintillation light and reflects light towards the central readout, enhancing
the detection efficiency. The detector holder plates are made out of PEN,
increasing the background suppression efficiency in the close vicinity of the
detectors. They scintillate when traversed by ionizing radiation and shift LAr
scintillation light to the visible blue.

To understand and model the response of the LAr detector of LEGEND-200 and
other experiments, the optical properties of the WLS must be known precisely.
This includes, among others, the emission spectrum and the wavelength shifting
efficiency. For LAr-based experiments, this means, that the materials need to
be characterized at $87\text{\,}\mathrm{K}$ and with VUV light excitation,
i.e., $128\text{\,}\mathrm{nm}$.

##  2 Cryogenic VUV spectrofluorometer

A dedicated experimental setup was designed and commissioned to characterize
WLS at temperatures ranging from
$300\text{\,}\mathrm{K}87\text{\,}\mathrm{K}$, with excitation wavelengths
down to $128\text{\,}\mathrm{nm}$. The VUV spectrofluorometer is schematically
depicted in Figure 1.

Figure 1: Schematic drawing of the VUV spectrofluorometer. The Roman numerals
indicate components referenced in the text.

In the sample chamber, a copper holder accommodates a $30\times
30\,\text{mm}^{2}$ WLS sample and is equipped with a $30\text{\,}\mathrm{W}$
polyimide heater on its back (indicated in Figure 1 as (i)). The holder,
positioned on the second stage of a two-stage Gifford McMahon coldhead, is
encased by a copper cryo-shield (ii), leaving openings for excitation light,
light detection, and cabling. A copper beaker filled with silica gel is
connected to a copper cold finger, which is partially submerged in a liquid
nitrogen dewar outside the chamber (iii). The temperature at the sample holder
recorded by Pt-100 sensors is the input for a PID controlling the heater’s
power (iv)[6].

A Hamamatsu L15094 deuterium lamp
($\lambda_{\text{Lamp}}=$120\text{\,}\mathrm{nm}$-$190\text{\,}\mathrm{nm}$$)
(v) and a Marktech Optoelectronics MTE310H33-UV LED
($\lambda_{\text{LED}}=$310\pm 10\text{\,}\mathrm{nm}$$) (vi) are used to
excite the sample. The deuterium lamp is mounted to an Acton VM-502 VUV-IR
monochromator (vii), which is attached to the sample chamber, enabling
excitation light at specific wavelengths ($20\text{\,}\mathrm{nm}$ width). The
deuterium lamp degrades in light output over time[7]. Therefore, a VUV-
sensitive McPherson photodiode (viii) is used to monitor the output of the
deuterium lamp. The LED is operated by a constant current LED driver and is
mounted on the cryo-shield.

For light collection, a collimating lens directs some of the shifted light to
an optical fiber, guiding it to an OceanOptics QE65000 UV-IR spectrometer
(ix)[7]. The lens is perpendicular to the excitation beam. By rotating the
sample holder, the angle of incidence and detection can be adjusted between
$15\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$75\text{\,}\mathrm{\SIUnitSymbolDegree}$ in
$15\text{\,}\mathrm{\SIUnitSymbolDegree}$ steps. With the setup, it is
possible to measure the relative wavelength shifting efficiency and emission
spectra at different temperatures or compare the properties between different
samples in the same configuration.

To prevent attenuation of VUV light on air, the entire system is evacuated to
a vacuum of approximately ${10}^{-7}\text{\,}\mathrm{mbar}$. When cooling the
sample, it is possible for residual humidity to form a thin ice layer on the
surface of the sample. This surface effect has been shown to significantly
reduce the amount of VUV light reaching the sample, affecting the measurement
[8]. To counter this effect, a two-step mitigation strategy is used: first,
the cryo-shield and beaker filled with silica gel are cooled before the sample
to adsorb the remaining humidity. Second, the UV LED is used as a
normalization for the measurements as it is unaffected by the ice layer.

##  3 Characterization of wavelength shifters

The WLSR sample used in this study is from the same material installed in
LEGEND-200 WLSR. The sample comprises of a $600\text{\,}\mathrm{nm}$ TPB layer
coated on $250\text{\,}\mathrm{nm}$ Tetratex®, backed by a copper foil.
Further details on how the WLSR was coated and details of the commissioning
can be found in [5].

The PEN sample was cut from the same material installed in LEGEND-200. The
material was formed using injection molding of commercially available PEN
pellets from Teĳin-DuPont (TN-8065 SC). Additional information on the
material’s production and characterization is presented in [4]. The amorphous
PEN sample guides its emission light to its edges, resulting in a low signal-
to-noise ratio in the VUV spectrofluorometer. To enhance the amount of
detected light, the surface of the PEN sample was sanded using a 600-grit
diamond sanding block, and a copper foil was placed as a reflector and thermal
coupling to the sample holder behind the sample.

For each measurement, the obtained spectrum is corrected for the known
response of the light detection system and the dark noise of the spectrometer.
Furthermore, the data is refined by excluding noisy pixels. The stability of
the deuterium lamp is monitored using the VUV-sensitive photodiode. Three
sources of uncertainties are considered in the data analysis: 1) system
instability, 2) uncertainty in the dark noise level, and 3) uncertainty in the
response of the light detection system. To obtain the relative WLSE, the WLSR
and PEN spectra are integrated between
$370\text{\,}\mathrm{nm}550\text{\,}\mathrm{nm}$ and
$370\text{\,}\mathrm{nm}650\text{\,}\mathrm{nm}$, respectively.

##  4 Results and discussion

(a) at $128\text{\,}\mathrm{nm}$ excitation.

(b) at $310\text{\,}\mathrm{nm}$ excitation.

Figure 2: Emission spectra of the WLSR sample for temperatures between
$293\text{\,}\mathrm{K}87\text{\,}\mathrm{K}$.

The obtained spectra for the characterization of the WLSR sample are shown in
Figure 2. The vibronic structures of the TPB molecule at
$405\text{\,}\mathrm{nm}450\text{\,}\mathrm{nm}$ emerge at lower temperatures.
In Table 1, the relative wavelength shifting efficiency (WLSE) of the
cryogenic measurements compared to the measurement at room temperature (RT) is
shown. We measured a strong increase in WLSE of $54\pm 4\text{\,}\mathrm{\char
37\relax}$ at LAr temperature ($87\text{\,}\mathrm{K}$) for excitation with
LAr scintillation light wavelength ($128\text{\,}\mathrm{nm}$). The intrinsic
quantum efficiency (QE) of TPB at RT for VUV excitation was measured by Benson
et al. to be $0.60\pm 0.04$[9]. Multiplied with the relative increase measured
in this work, one gets a QE of approximately $0.92\pm 0.08$, which is
compatible with the QE of the LEGEND-200 WLSR obtained by Araujo et al. of
$0.85\pm 0.05$(stat)$\pm 0.06$(syst.)[10].

Table 1: Relative WLSE of the WLSR and PEN sample for cryogenic temperatures. Temperature [K] | Relative WLSE compared to RT [%]
---|---
WLSR | PEN
$\lambda_{\text{exc}}=$128\text{\,}\mathrm{nm}$$ | $\lambda_{\text{exc}}=$310\text{\,}\mathrm{nm}$$ | $\lambda_{\text{exc}}=$128\text{\,}\mathrm{nm}$$ | $\lambda_{\text{exc}}=$310\text{\,}\mathrm{nm}$$
250 | $112\pm 4$ | $106\pm 2$ | $112\pm 2$ | $110\pm 2$
200 | $118\pm 4$ | $113\pm 2$ | $123\pm 3$ | $122\pm 2$
150 | $130\pm 4$ | $127\pm 2$ | $134\pm 4$ | $137\pm 3$
100 | $147\pm 5$ | $142\pm 2$ | $134\pm 4$ | $152\pm 2$
87 | $154\pm 5$ | $144\pm 2$ | $137\pm 4$ | $152\pm 3$

Our results are compared with previous studies of TPB. Ellingwood et al.
investigated the emission of TPB excited by UV LED with a wavelength of
$\lambda_{\text{exc}}=$285\text{\,}\mathrm{nm}$$. For their experiment, a
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick TPB layer was coated on a
$5\text{\,}\mathrm{mm}$ thick PMMA substrate[11]. Francini et al.
characterized a TPB coating of approximately
$1.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thickness on polymeric reflector
(Vikuiti™) deposited by vacuum evaporation [12]. The emission spectrum at
$87\text{\,}\mathrm{K}$ is compared to these studies in 3(a).

(a) Emission Spectra at $87\text{\,}\mathrm{K}$.

(b) Relative WLSE compared to RT for temperatures of
$300\text{\,}\mathrm{K}87\text{\,}\mathrm{K}$.

Figure 3: Comparison to the results obtained by Ellingwood et al.[11], and
Francini et al.[12].

The spectrum we detected is in good agreement with the results from [11] and
[12]. The minor differences in the spectra might originate in differences in
the TPB vacuum evaporation procedure[7], in differences in the
reflectance/transmittance of the backing material, or the lack of calibrated
response of the detection system used in [11] and [12].

In 3(b), the relative WLSE between RT and $87\text{\,}\mathrm{K}$ is compared.
The relative increase we measured both with $128\text{\,}\mathrm{nm}$ and
$310\text{\,}\mathrm{nm}$ excitation exceeds the increase reported by [11] and
[12] by roughly a factor of two. The measurements by Francini et al. might
suffer strongly from the above-mentioned ice layer formation, lowering the
increase in WLSE[8]. Alternatively, the quality of the TPB layer and the type
of substrate may influence the performance at cold temperatures differently.

(a) at $128\text{\,}\mathrm{nm}$ excitation.

(b) at $310\text{\,}\mathrm{nm}$ excitation.

Figure 4: Emission spectra of the PEN sample for temperatures between
$293\text{\,}\mathrm{K}87\text{\,}\mathrm{K}$.

The emission spectra of the PEN sample for temperatures between
$300\text{\,}\mathrm{K}$ and $87\text{\,}\mathrm{K}$ are displayed in Figure
4. A slight shift in the peak wavelength from $430\text{\,}\mathrm{nm}$ for
$293\text{\,}\mathrm{K}$ to $425\text{\,}\mathrm{nm}$ for
$87\text{\,}\mathrm{K}$ can be seen. Additionally, the monomeric emission at
$575\text{\,}\mathrm{nm}$ is minimally enhanced at cold temperatures. Both
effects are in agreement with previous studies of PEN using
$300\text{\,}\mathrm{nm}$ excitation[13].

The relative increase in WLSE compared to $293\text{\,}\mathrm{K}$ is
presented in Table 1. One can see a good agreement between the increase for
both excitation wavelengths down to $150\text{\,}\mathrm{K}$. Below, the WLSE
increases significantly less for the VUV excitation. Because of its larger
mass compared to the WLSR sample and its poor thermal conductivity, the PEN
sample needs more time to cool down to $87\text{\,}\mathrm{K}$. Therefore, the
surface effect, which becomes significant below $150\text{\,}\mathrm{K}$ for
${10}^{-7}\text{\,}\mathrm{mbar}$[14], impacts the measurement with
$128\text{\,}\mathrm{nm}$ much stronger than for the WLSR measurement. Hence,
one can assume that the WLSE of PEN at LAr temperature
($87\text{\,}\mathrm{K}$) relative to $293\text{\,}\mathrm{K}$ for excitation
with LAr scintillation light wavelength ($128\text{\,}\mathrm{nm}$) is equal
or higher than the measured value for $310\text{\,}\mathrm{nm}$ excitation of
$149\pm 1\text{\,}\mathrm{\char 37\relax}$.

##  5 Conclusion

We present the necessity and challenge of the characterization of WLS for the
application in LAr-based experiments like LEGEND-200. For the
characterization, we commissioned a VUV spectrofluorometer capable of
measuring the emission spectrum and relative WLSE of wavelength shifting
samples for excitation with LAr scintillation wavelength
($128\text{\,}\mathrm{nm}$) at temperatures from RT down to LAr temperature
($87\text{\,}\mathrm{K}$). For the TPB-based wavelength shifting reflector
featured in LEGEND-200, the WLSE at $87\text{\,}\mathrm{K}$ was found to be
increased by $54\pm 5\text{\,}\mathrm{\char 37\relax}$ compared to RT. The
WLSE of PEN, used in LEGEND-200 for the optically active detector holders,
increased by at least $49\pm 1\text{\,}\mathrm{\char 37\relax}$ by cooling
from $293\text{\,}\mathrm{K}87\text{\,}\mathrm{K}$. To the knowledge of the
authors, this work is the first characterization of the emission spectrum and
relative WLSE of amorphous PEN at cryogenic temperatures with VUV and UV
excitation.

## Acknowledgments

This research is supported by the DFG through the Excellence Cluster ORIGINS
and the SFB1258.

## References

  * [1] P.A. Amaudruz, M. Baldwin, M. Batygov, B. Beltran, C.E. Bina, D. Bishop et al., _Design and construction of the DEAP-3600 dark matter detector_ , _Astroparticle Physics_ 108 (2019) 1.
  * [2] C.E. Aalseth, F. Acerbi, P. Agnes, I.F.M. Albuquerque, T. Alexander, A. Alici et al., _DarkSide-20k: A 20 tonne two-phase LAr TPC for direct dark matter detection at LNGS_ , _The European Physical Journal Plus_ 133 (2018) 131.
  * [3] N. Abgrall, I. Abt, M. Agostini, A. Alexander, C. Andreoiu, G.R. Araujo et al., _LEGEND-1000 Preconceptual Design Report_ , July, 2021.  10.48550/arXiv.2107.11462.
  * [4] Y. Efremenko, M. Febbraro, F. Fischer, M.G. Corominas, K. Gusev, B. Hackett et al., _Production and validation of scintillating structural components from low-background Poly(ethylene naphthalate)_ , _Journal of Instrumentation_ 17 (2022) P01010.
  * [5] M. Schwarz, P. Krause, A. Leonhardt, L. Papp, S. Schönert, C. Wiesinger et al., _Liquid Argon Instrumentation and Monitoring in LEGEND-200_ , _EPJ Web of Conferences_ 253 (2021) 11014.
  * [6] M. Goldbrunner, _Characterization of PEN-based Wavelength Shifting Reflectors for the LEGEND Experiment_ , Feb., 2023.  10.13140/RG.2.2.26723.53288.
  * [7] A. Leonhardt, _Characterization of Wavelength Shifters for Rare-Event Search Experiments with a VUV Spectrofluorometer_ , July, 2021.  10.13140/RG.2.2.25067.87849.
  * [8] A. Neumeier, T. Dandl, A. Himpsl, M. Hofmann, L. Oberauer, W. Potzel et al., _Attenuation measurements of vacuum ultraviolet light in liquid argon revisited_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 800 (2015) 70.
  * [9] C. Benson, G.D. Orebi Gann and V. Gehman, _Measurements of the intrinsic quantum efficiency and absorption length of tetraphenyl butadiene thin films in the vacuum ultraviolet regime_ , _The European Physical Journal C_ 78 (2018) 329.
  * [10] G.R. Araujo, L. Baudis, N. McFadden, P. Krause, S. Schönert and V.H.S. Wu, _R &D of wavelength-shifting reflectors and characterization of the quantum efficiency of tetraphenyl butadiene and polyethylene naphthalate in liquid argon_, _The European Physical Journal C_ 82 (2022) 442.
  * [11] E. Ellingwood, H. Benmansour, Q. Hars, J. Hucker, V. Pereymak, J. Corning et al., _Ultraviolet-induced fluorescence of poly(methyl methacrylate) compared to 1,1,4,4-tetraphenyl-1,3-butadiene down to 4 k_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 1039 (2022) 167119.
  * [12] R. Francini, R.M. Montereali, E. Nichelatti, M.A. Vincenti, N. Canci, E. Segreto et al., _Vuv-vis optical characterization of tetraphenyl-butadiene films on glass and specular reflector substrates from room to liquid argon temperature_ , _Journal of Instrumentation_ 8 (2013) P09006.
  * [13] D. Mary, M. Albertini and C. Laurent, _Understanding optical emissions from electrically stressed insulating polymers: electroluminescence in poly(ethylene terephthalate) and poly(ethylene 2,6-naphthalate) films_ , _Journal of Physics D: Applied Physics_ 30 (1997) 171.
  * [14] J.M. Labello, _Water Ice Films in Cryogenic Vacuum Chambers_ , _Doctoral Dissertations_ (2011) .