---
abstract: |
    In this investigation we study extreme vortex states defined as incompressible velocity fields with prescribed enstrophy $\E_0$ which maximize the instantaneous rate of growth of enstrophy $d\E/dt$. We provide [an analytic]{} characterization of these extreme vortex states in the limit of vanishing enstrophy $\E_0$ and, in particular, show that the Taylor-Green vortex is in fact a local maximizer of $d\E / dt$ [in this limit]{}. For finite values of enstrophy, the extreme vortex states are computed numerically by solving a constrained variational optimization problem using a suitable gradient method. In combination with a continuation approach, this allows us to construct an entire family of maximizing vortex states parameterized by their enstrophy. We also confirm the findings of the seminal study by [@ld08] that these extreme vortex states saturate (up to a numerical prefactor) the fundamental bound $d\E / dt < C \, \E^3$, for some constant $C >
      0$. The time evolution corresponding to these extreme vortex states leads to a larger growth of enstrophy than the growth achieved by any of the commonly used initial conditions with the same enstrophy $\E_0$. However, based on several different diagnostics, there is no evidence of any tendency towards singularity formation in finite time. Finally, we discuss possible physical reasons why the initially large growth of enstrophy is not sustained for longer times.
author:
- |
    Diego Ayala$^{1,2}$ and Bartosz Protas$^{2,}$[^1]\
    \
    $^1$ Department of Mathematics, University of Michigan,\
    Ann Arbor, MI 48109, USA\
    \
    $^2$ Department of Mathematics and Statistics, McMaster University\
    Hamilton, Ontario, L8S 4K1, Canada
title: Extreme Vortex States and the Growth of Enstrophy in 3D Incompressible Flows
---

Keywords: Navier-Stokes equations; Extreme Behaviour; Variational methods; Vortex Flows

Introduction {#sec:intro}
============

The objective of this investigation is to study three-dimensional (3D) flows of viscous incompressible fluids which are constructed to exhibit extreme growth of enstrophy. It is motivated by the question whether the solutions to the 3D incompressible Navier-Stokes system on unbounded or periodic domains corresponding to smooth initial data may develop a singularity in finite time [@d09]. By formation of a “singularity” we mean the situation when some norms of the solution corresponding to smooth initial data have become unbounded after a finite time. This so-called “blow-up problem” is one of the key open questions in mathematical fluid mechanics and, in fact, its importance for mathematics in general has been recognized by the Clay Mathematics Institute as one of its “millennium problems” [@f00]. Questions concerning global-in-time existence of smooth solutions remain open also for a number of other flow models including the 3D Euler equations [@gbk08] and some of the “active scalar” equations [@k10].

While the blow-up problem is fundamentally a question in mathematical analysis, a lot of computational studies have been carried out since the mid-’90s in order to shed light on the hydrodynamic mechanisms which might lead to singularity formation in finite time. Given that such flows evolving near the edge of regularity involve formation of very small flow structures, these computations typically require the use of state-of-the-art computational resources available at a given time. The computational studies focused on the possibility of finite-time blow-up in the 3D Navier-Stokes and/or Euler system include [@bmonmu83; @ps90; @b91; @k93; @p01; @bk08; @oc08; @o08; @ghdg08; @gbk08; @h09; @opc12; @bb12; @opmc14], all of which considered problems defined on domains periodic in all three dimensions. Recent investigations by [@dggkpv13; @k13; @gdgkpv14; @k13b] focused on the time evolution of vorticity moments and compared it with the predictions derived from analysis based on rigorous bounds. We also mention the studies by [@mbf08] and [@sc09], along with the references found therein, in which various complexified forms of the Euler equation were investigated. The idea of this approach is that, since the solutions to complexified equations have singularities in the complex plane, singularity formation in the real-valued problem is manifested by the collapse of the complex-plane singularities onto the real axis.

Overall, the outcome of these investigations is rather inconclusive: while for the Navier-Stokes flows most of recent computations do not offer support for finite-time blow-up, the evidence appears split in the case of the Euler system. In particular, the recent studies by [@bb12] and [@opc12] hinted at the possibility of singularity formation in finite time. In this connection we also mention the recent investigations by [@lh14a; @lh14b] in which blow-up was observed in axisymmetric Euler flows in a bounded (tubular) domain.

A common feature of all of the aforementioned investigations was that the initial data for the Navier-Stokes or Euler system was chosen in an ad-hoc manner, based on some heuristic arguments. On the other hand, in the present study we pursue a fundamentally different approach, proposed originally by [@ld08] and employed also by [@ap11a; @ap13a; @ap13b] for a range of related problems, in which the initial data leading to the most singular behaviour is sought systematically via solution of a suitable variational optimization problem. We carefully analyze the time evolution induced by the extreme vortex states first identified by [@ld08] and compare it to the time evolution corresponding to a number of other candidate initial conditions considered in the literature [@bmonmu83; @k93; @p01; @cb05; @opc12]. We demonstrate that the Taylor-Green vortex, studied in the context of the blow-up problem by [@tg37; @bmonmu83; @b91; @cb05], is in fact a particular member of the family of extreme vortex states maximizing the instantaneous rate of enstrophy production in the limit of vanishing enstrophy. In addition, based on these findings, we identify the set of initial data, parameterized by its energy and enstrophy, for which one can a priori guarantee global-in-time existence of smooth solutions. This result therefore offers a physically appealing interpretation of an “abstract” mathematical theorem concerning global existence of classical solutions corresponding to “small” initial data [@lady69]. We also emphasize that, in order to establish a direct link with the results of the mathematical analysis discussed below, in our investigation we follow a rather different strategy than in most of the studies referenced above. While these earlier studies relied on data from a relatively small number of simulations performed at a high (at the given time) resolution, in the present investigation we explore a broad range of cases, each of which is however computed at a more moderate resolution (or, equivalently, Reynolds number). With such an approach to the use of available computational resources, we are able to reveal trends resulting from the variation of parameters which otherwise would be hard to detect. Systematic computations conducted in this way thus allow us to probe the sharpness of the mathematical analysis relevant to the problem.

The question of regularity of solution to the Navier-Stokes system is [usually addressed using “energy” methods which rely on finding upper bounds (with respect to time) on certain quantities of interest, typically taken as suitable Sobolev norms of the solution. A key intermediate step is obtaining bounds on the rate of growth of the quantity of interest, a problem which can be studied with ODE methods.]{} While [for the Navier-Stokes system]{} different norms of the velocity gradient or vorticity can be used to study the regularity of solutions, the use of enstrophy $\E$ (see equation ) is privileged by the well-known result of @ft89, where it was established that if the uniform bound $$\label{eq:RegCrit_FoiasTemam}
\mathop{\sup}_{0 \leq t \leq T} \E(\u(t))  < \infty$$ holds, then the regularity of the solution $\u(t)$ is guaranteed up to time $T$ (to be precise, the solution remains in a suitable Gevrey class). From the computational point of view, the enstrophy $\E(t) := \E(\u(t))$ is thus a convenient indicator of the regularity of solutions, because in the light of , singularity formation must manifest itself by the enstrophy becoming infinite.

While characterization of the maximum possible finite-time growth of enstrophy in the 3D Navier-Stokes flows is the ultimate objective of this research program, analogous questions can also be posed in the context of more tractable problems involving the one-dimensional (1D) Burgers equation and the two-dimensional (2D) Navier-Stokes equation. Although global-in-time existence of the classical (smooth) solutions is well known for both these problems [@kl04], questions concerning the sharpness of the corresponding estimates for the instantaneous and finite-time growth of various quantities are relevant, because these estimates are obtained using essentially the same methods as [employed to derive]{} their 3D counterparts. Since in 2D flows on unbounded or periodic domains the enstrophy may not increase ($d\E/dt \leq 0$), the relevant quantity in this case is the palinstrophy $\P({\mathbf{u}}) := \frac{1}{2}\int_\Omega |
\bnabla{\bm{\omega}}({\mathbf{x}},t) |^2 \,d{\mathbf{x}}$, where ${\bm{\omega}}:={\bnabla\times}{\mathbf{u}}$ is the vorticity (which reduces to a pseudo-scalar in 2D). Different questions concerning sharpness of estimates addressed in our research program are summarized together with the results obtained to date in Table \[tab:estimates\]. We remark that the best finite-time estimate for the 1D Burgers equation was found [*not*]{} to be sharp using the initial data obtained from both the instantaneous and the finite-time variational optimization problems [@ap11a]. On the other hand, in 2D the bounds on both the *instantaneous* and *finite-time* growth of palinstrophy were found to be sharp and, somewhat surprisingly, both estimates were realized by the same family of incompressible vector fields parameterized by energy $\K$ and palinstrophy $\P$, obtained as the solution of an [*instantaneous*]{} optimization problem [@ap13a]. It is worth mentioning that while the estimate for the instantaneous rate of growth of palinstrophy $d\P/dt \leq C\K^{1/2}\P^{3/2}/\nu$ (see Table \[tab:estimates\]) was found to be sharp with respect to variations in palinstrophy, the estimate is in fact not sharp with respect to the prefactor $C_{{\mathbf{u}},\nu} = \K^{1/2}/\nu$ [@ados16], with the correct prefactor being of the form $\widetilde{C}_{{\mathbf{u}},\nu} =
\sqrt{\log\left(\K^{1/2}/\nu\right)}$. We add that what distinguishes the 2D problem, in regard to both the instantaneous and finite-time bounds, is that the RHS of these bounds are expressed in terms of two quantities, namely, energy $\K$ and enstrophy $\E$, in contrast to the enstrophy alone appearing in the 1D and 3D estimates. As a result, the 2D instantaneous optimization problem had to be solved subject to [ *two*]{} constraints.

In the present investigation we advance the research program summarized in Table \[tab:estimates\] by assessing to what extent the finite-time growth of enstrophy predicted by the analytic estimates and can be actually realized by flow evolution starting from different initial conditions, including the extreme vortex states found by [@ld08] to saturate the instantaneous estimate . The key finding is that, at least for the range of modest enstrophy values we considered, the growth of enstrophy corresponding to this initial data, which has the form of two colliding axisymmetric vortex rings, is rapidly depleted and there is no indication of singularity formation in finite time. Thus, should finite-time singularity be possible in the Navier-Stokes system, it is unlikely to result from initial conditions instantaneously maximizing the rate of growth of enstrophy. We also provide a comprehensive characterization of the extreme vortex states which realize estimate together with the resulting flow evolutions.

The structure of the paper is as follows: in the next section we present analytic estimates [on]{} the instantaneous and finite-time growth of enstrophy in 3D flows. In §\[sec:3D\_InstOpt\] we formulate the variational optimization problems which will be solved to find the vortex states with the largest rate of enstrophy production and in §\[sec:3D\_InstOpt\_E0to0\] we provide an asymptotic representation for these optimal states in the limit of vanishing enstrophy. In §\[sec:3D\_InstOpt\_E\] we present numerically computed extreme vortex states corresponding to intermediate and large enstrophy values, while in §\[sec:timeEvolution\] we analyze the temporal evolution corresponding to different initial data in order to compare it with the predictions of estimates and . Our findings are discussed in §\[sec:discuss\], whereas conclusions and outlook are deferred to §\[sec:final\].

[l|c|c]{} &

[Estimate]{}\

&

[Realizability ]{}

\
& &\
& &\
& $\frac{d\P(t)}{dt}  \le -\nu\frac{\P^2}{\E} + \frac{C_1}{\nu} \E\,\P$\
$\frac{d\P(t)}{dt} \le \frac{C_2}{\nu} \K^{1/2}\P^{3/2}$ &\

2D Navier-Stokes\
finite-time

& $\max_{t>0} \P(t) \le \P_0 + \frac{C_1}{2\nu^2}\E_0^2$\
$\max_{t>0} \P(t) \le \left(\P_0^{1/2} + \frac{C_2}{4\nu^2}\K_0^{1/2}\E_0\right)^2$ &\
& &\

3D Navier-Stokes\
finite-time

& $\E(t) \le \frac{\E(0)}{\sqrt{1 - 4 \frac{C \E(0)^2}{\nu^3} t}}$\
$\frac{1}{\E(0)} - \frac{1}{\E(t)} \leq  \frac{27}{(2\pi\nu)^4}\left[\K(0) - \K(t) \right]$ &\

Bounds on the Growth of Enstrophy in 3D Navier-Stokes Flows {#sec:Bounds3DNS}
===========================================================

We consider the incompressible Navier-Stokes system defined on the 3D unit cube $\Omega = [0,1]^3$ with periodic boundary conditions

\[eq:NSE3D\] $$\begin{aligned}
{2}
\partial_t{\mathbf{u}}+ {\mathbf{u}}\cdot\bnabla{\mathbf{u}}+ \bnabla p - \nu{\Delta}{\mathbf{u}}& = 0 & &\qquad\mbox{in} \,\,\Omega\times(0,T), \\
\bnabla\cdot{\mathbf{u}}& = 0 & & \qquad\mbox{in} \,\,\Omega\times[0,T), \\
{\mathbf{u}}({\mathbf{x}},0) & = {\mathbf{u}}_0({\mathbf{x}}), &   &\end{aligned}$$

where the vector ${\mathbf{u}}= [u_1, u_2, u_3]$ is the velocity field, $p$ is the pressure and $\nu>0$ is the coefficient of kinematic viscosity (hereafter we will set $\nu=0.01$ which is the same value as used in the seminal study by [@ld08]). The velocity gradient $\bnabla{\mathbf{u}}$ is the tensor with components $[\bnabla{\mathbf{u}}]_{ij} =
  \partial_j u_i$, $i,j=1,2,3$. The fluid density $\rho$ is assumed to be constant and equal to unity ($\rho=1$). The relevant properties of solutions to system can be studied using energy methods, with the energy $\K({\mathbf{u}})$ and its rate of growth given by $$\begin{aligned}
  \K({\mathbf{u}}) & := & \frac{1}{2}\int_\Omega |{\mathbf{u}}({\mathbf{x}},t)|^2 \,d{\mathbf{x}}, \label{eq:EnerDef_3D}\\
  \frac{d\K({\mathbf{u}})}{dt} & = & -\nu\int_\Omega |\nabla{\mathbf{u}}|^2 \, d{\mathbf{x}}, \label{eq:dK/dt_3D}\end{aligned}$$ where “$:=$” means “equal to by definition”. The enstrophy $\E({\mathbf{u}})$ and its rate of growth are given by $$\begin{aligned}
\E({\mathbf{u}}) & := & \frac{1}{2}\int_\Omega | {\bnabla\times}{\mathbf{u}}({\mathbf{x}},t) |^2 \,d{\mathbf{x}}, \label{eq:EnsDef_3D}\\
\frac{d\E({\mathbf{u}})}{dt} & = & -\nu\int_\Omega |{\Delta}{\mathbf{u}}|^2\,d{\mathbf{x}}+ 
\int_{\Omega} {\mathbf{u}}\cdot\nabla{\mathbf{u}}\cdot{\Delta}{\mathbf{u}}\, d{\mathbf{x}}=: \R({\mathbf{u}}). \label{eq:dEdt}\end{aligned}$$ For incompressible flows with periodic boundary conditions we also have the following identity [@dg95] $$\int_{\Omega} |{\bnabla\times}{\mathbf{u}}|^2\,d{\mathbf{x}}= \int_{\Omega} |\nabla{\mathbf{u}}|^2\,d{\mathbf{x}}.
\label{eq:duL2}$$ Hence, combining –, the energy and enstrophy satisfy the system of ordinary differential equations

$$\begin{aligned}
\frac{d\K({\mathbf{u}})}{dt} & =  -2\nu\E({\mathbf{u}}), \label{eq:dKdt_system}\\
\frac{d\E({\mathbf{u}})}{dt} & =  \R({\mathbf{u}}). \label{eq:dEdt_system}\end{aligned}$$

A standard approach at this point is to try to upper-bound $d\E / dt$ and using standard techniques of functional analysis it is possible to obtain the following well-known estimate in terms of $\K$ and $\E$ [@d09] $$\label{eq:dEdt_estimate_KE}
\frac{d\E}{dt} \leq -\nu \frac{\E^2}{\K} + \frac{c}{\nu^3} \E^3$$ for $c$ an absolute constant. A related estimate expressed entirely in terms of the enstrophy $\E$ is given by $$\frac{d\E}{dt} \leq \frac{27}{8\,\pi^4\,\nu^3} \E^3. 
\label{eq:dEdt_estimate_E}$$ By simply integrating the differential inequality in with respect to time we obtain the finite-time bound $$\E(t) \leq \frac{\E(0)}{\sqrt{1 - \frac{27}{4\,\pi^4\,\nu^3}\,\E(0)^2\, t}}
\label{eq:Et_estimate_E0}$$ which clearly becomes infinite at time $t_0 = 4\,\pi^4\,\nu^3 /
[27\,\E(0)^2]$. Thus, based on estimate , it is not possible to establish the boundedness of the enstrophy $\E(t)$ globally in time and hence the regularity of solutions. Therefore, the question about the finite-time singularity formation can be recast in terms of whether or not estimate can be saturated. By this we mean the existence of initial data with enstrophy $\E_0 := \E(0)> 0$ such that the resulting time evolution realizes the largest growth of enstrophy $\E(t)$ allowed by the right-hand side (RHS) of estimate . A systematic search for such most singular initial data using variational optimization methods is the key theme of this study. Although different notions of sharpness of an estimate can be defined, e.g., sharpness with respect to constants or exponents in the case of estimates in the form of power laws, the precise notion of sharpness considered in this study is the following

\[def:NotionSharpness\] Given a parameter $p\in\mathbb{R}$ and maps $f,g:\mathbb{R}\to\mathbb{R}$, the estimate $$f(p) \leq g(p)$$ is declared sharp in the limit $p \to p_0\in\mathbb{R}$ if and only if $$\lim_{p \to p_0} \frac{f(p)}{g(p)} \sim \beta, \quad \beta \in \RR.$$

From this definition, the sharpness of estimates in the form $g(p) =
C\, p^{\alpha}$ for some $C \in \RR_+$ and $\alpha \in \RR$ can be addressed in the limit $p \rightarrow \infty$ by studying the adequacy of the exponent $\alpha$.

The question of sharpness of estimate was addressed in the seminal study by [@ld08], see also [@l06], who constructed a family of divergence-free velocity fields saturating this estimate. More precisely, these vector fields were parameterized by their enstrophy and for sufficiently large values of $\E$ the corresponding rate of growth $d\E/dt$ was found to be proportional to $\E^3$. Therefore, in agreement with definition \[def:NotionSharpness\], estimate was declared sharp up to a numerical prefactor. However, the sharpness of the instantaneous estimate alone does not allow us to conclude about the possibility of singularity formation, because for this situation to occur, a sufficiently large enstrophy growth rate would need to be sustained over a [*finite*]{} time window $[0,t_0)$. In fact, assuming the instantaneous rate of growth of enstrophy in the form $d\E / dt =
C \, \E^{\alpha}$ for some $C>0$, any exponent $\alpha > 2$ will produce blow-up of $\E(t)$ in finite time if the rate of growth is sustained. The fact that there is no blow-up for $\alpha \le
  2$ follows from Grönwall’s lemma and the fact that one factor of $\E$ in can be bounded in terms of the initial energy using as follows $$\int_0^t \E(s)\, ds = \frac{1}{2\nu} \left[ \K(0) - \K(t)\right] \leq \frac{1}{2\nu} \K(0).
\label{eq:Kt}$$ This relation also leads to an alternative form of the estimate for the finite-time growth of enstrophy, namely $$\begin{aligned}
\frac{d\E}{dt} & \leq  \frac{27}{8\pi^4\nu^3}\E^3\quad\Longrightarrow \nonumber \\
\int_{\E(0)}^{\E(t)}\E^{-2}\,d\E & \leq  \frac{27}{8\pi^4\nu^3}\int_0^t\E(s)\,ds \quad\Longrightarrow \nonumber \\
\frac{1}{\E(0)} - \frac{1}{\E(t)} & \leq  \frac{27}{(2\pi\nu)^4}\left[\K(0) - \K(t) \right] \label{eq:Evs_t_fixE}\end{aligned}$$ which is more convenient than from the computational point of view and will be used in the present study. We note, however, that since the RHS of this inequality cannot be expressed entirely in terms of properties of the initial data, this is [*not*]{} in fact an a priori estimate. Estimate also allows us to obtain a condition on the size of the initial data, given in terms of its energy $\K(0)$ and enstrophy $\E(0)$, which guarantees that smooth solutions will exist globally in time, namely, $$\label{eq:Cond_for_globalReg}
\mathop{\max}_{t \geq 0} \E(t) \leq \frac{\E(0)}{1 - \frac{27}{(2\pi\nu)^4}\K(0)\E(0)}$$ from which it follows that $$\label{eq:K0E0}
\K(0)\E(0) < \frac{(2\pi\nu)^4}{27}.$$ Thus, flows with energy and enstrophy satisfying inequality are guaranteed to be smooth for all time, in agreement with [the]{} regularity results [available under the assumption of]{} small initial data [[@lady69]]{}.

Instantaneously Optimal Growth of Enstrophy {#sec:3D_InstOpt}
===========================================

Sharpness of instantaneous estimate , in the sense of definition \[def:NotionSharpness\], can be probed by constructing a family of “extreme vortex states” ${\widetilde{\mathbf{u}}_{\E_0}}$ which, for each $\E_0 > 0$, have prescribed enstrophy $\E({\widetilde{\mathbf{u}}_{\E_0}}) =
  \E_0$ and produce the largest possible rate of growth of enstrophy $\R({\widetilde{\mathbf{u}}_{\E_0}})$. Given the form of , the fields ${\widetilde{\mathbf{u}}_{\E_0}}$ can be expected to exhibit (at least piecewise) smooth dependence on $\E_0$ and we will refer to the mapping $\E_0
  \longmapsto {\widetilde{\mathbf{u}}_{\E_0}}$ as a “[maximizing]{} branch”. Thus, information about the sharpness of estimate can be deduced by analyzing the relation $\E_0$ versus $\R({\widetilde{\mathbf{u}}_{\E_0}})$ obtained for a possibly broad range of enstrophy values. A [maximizing]{} branch is constructed by finding, for different values of $\E_0$, the extreme vortex states ${\widetilde{\mathbf{u}}_{\E_0}}$ as solutions of a variational optimization problem defined below.

Hereafter, $H^2(\Omega)$ will denote the Sobolev space of functions with square-integrable second derivatives endowed with the inner product [@af05] $$\forall\,\mathbf{z}_1, \mathbf{z}_2 \in H^2(\Omega) \qquad 
\Big\langle \mathbf{z}_1, \mathbf{z}_2 \Big\rangle_{H^2(\Omega)}
= \int_{\Omega} \mathbf{z}_1 \cdot \mathbf{z}_2 
+ \ell_1^2 \,\bnabla \mathbf{z}_1 \colon \bnabla \mathbf{z}_2
+ \ell_2^4 \,\Delta \mathbf{z}_1 \cdot \Delta \mathbf{z}_2  \, d{\mathbf{x}},  \label{eq:ipH2}$$ where $\ell_1,\ell_2\in \RR_+$ are parameters with the meaning of length scales (the reasons for introducing these parameters in the definition of the inner product will become clear below). The inner product in the space $L_2(\Omega)$ is obtained from by setting $\ell_1 = \ell_2 = 0$. The notation $H^2_0(\Omega)$ will refer to the Sobolev space $H^2(\Omega)$ of functions with zero mean. For every fixed value $\E_0$ of enstrophy we will look for a divergence-free vector field ${\widetilde{\mathbf{u}}_{\E_0}}$ maximizing the objective function $\R \; : \; H^2_0(\Omega) \rightarrow \RR$ defined in . We thus have the following

\[pb:maxdEdt\_E\] Given $\E_0\in\mathbb{R}_+$ and the objective functional $\R$ from equation , find $$\begin{aligned}
{\widetilde{\mathbf{u}}_{\E_0}}& = & \mathop{\arg\max}_{{\mathbf{u}}\in{\mathcal{S}_{\E_0}}} \, \R({\mathbf{u}}) \\
{\mathcal{S}_{\E_0}} & = & \left\{{\mathbf{u}}\in H_0^2(\Omega)\,\colon\,\nabla\cdot{\mathbf{u}}= 0, \; \E({\mathbf{u}}) = \E_0 \right\}\end{aligned}$$

which will be solved for enstrophy $\E_0$ spanning a broad range of values. This approach was originally proposed and investigated by [@ld08]. In the present study we extend and generalize these results by first showing how other fields considered in the context of the blow-up problem for both the Euler and Navier-Stokes system, namely the Taylor-Green vortex, also arise from variational problem \[pb:maxdEdt\_E\]. We then thoroughly analyze the time evolution corresponding to our extreme vortex states and compare it with the predictions of the finite-time estimates and . As discussed at the end of this section, some important aspects of our approach to solving problem \[pb:maxdEdt\_E\] are also quite different from the method adopted by [@ld08].

The smoothness requirement in the statement of problem \[pb:maxdEdt\_E\] (${\mathbf{u}}\in H_0^2(\Omega)$) follows from the definition of the objective functional $\R$ in equation , where both the viscous term $\nu\int_\Omega
|{\Delta}{\mathbf{u}}|^2\,d{\mathbf{x}}$ and the cubic term $\int_{\Omega}
{\mathbf{u}}\cdot\bnabla{\mathbf{u}}\cdot{\Delta}{\mathbf{u}}\, d{\mathbf{x}}$ contain derivatives of order up to two. The constraint manifold ${\mathcal{S}_{\E_0}}$ can be interpreted as an intersection of the manifold (a subspace) ${\mathcal{S}_{0}}
\in H_0^2(\Omega)$ of divergence-free fields and the manifold $\mathcal{S}'_{\E_0} \in H_0^2(\Omega)$ of fields with prescribed enstrophy $\E_0$. The structure of these constraint manifolds is reflected in the definition of the corresponding projections $\mathbb{P}_{\mathcal{S}}:H_0^2\to\mathcal{S}$ (without a subscript, $\mathcal{S}$ refers to a generic manifold) which is given for each of the two constraints as follows:

-   ([*div-free*]{})-constraint: the projection of a field ${\mathbf{u}}$ onto the subspace of solenoidal fields $\mathcal{S}_0$ is performed using the Helmholtz decomposition; accordingly, every zero-mean vector field ${\mathbf{u}}\in H_0^2(\Omega)$ can be decomposed uniquely as $${\mathbf{u}}= \bnabla\phi + {\bnabla\times}{\mathbf{A}},$$ where $\phi$ and ${\mathbf{A}}$ are scalar and vector potentials, respectively; it follows from the identity $\bnabla\cdot({\bnabla\times}{\mathbf{A}})\equiv0$, valid for any sufficiently smooth vector field ${\mathbf{A}}$, that the projection $\mathbb{P}_{\mathcal{S}_0}({\mathbf{u}})$ is given simply by ${\bnabla\times}{\mathbf{A}}$ and is therefore calculated as $$\mathbb{P}_{\mathcal{S}_0}({\mathbf{u}}) = {\mathbf{u}}- \bnabla\left[{\Delta}^{-1}(\bnabla\cdot{\mathbf{u}})\right], 
    \label{eq:Phodge}$$ where ${\Delta}^{-1}$ is the inverse Laplacian associated with the periodic boundary conditions; the operator $\mathbb{P}_{\mathcal{S}_0}$ is also known as the Leray-Helmholtz projector.

-   $(\E_0)$-constraint: the projection onto the manifold $\mathcal{S}'_{\E_0}$ is calculated by the normalization $$\label{eq:FixE0_3D}
    \mathbb{P}_{\mathcal{S}'_{\E_0}}({\mathbf{u}}) = \sqrt{\frac{\E_0}{\E\left({\mathbf{u}}\right)}}\,{\mathbf{u}}.$$

Thus, composing with , the projection onto the manifold ${\mathcal{S}_{\E_0}}$ defined in problem \[pb:maxdEdt\_E\] is constructed as $$\mathbb{P}_{\mathcal{S}_{\E_0}}({\mathbf{u}}) = \mathbb{P}_{\mathcal{S}'_{\E_0}}\Big(  \mathbb{P}_{\mathcal{S}_0} ({\mathbf{u}})\Big).
\label{eq:P}$$ This approach, which was already successfully employed by [@ap11a; @ap13a], allows one to enforce the enstrophy constraint essentially with the machine precision.

For a given value of $\E_0$, the maximizer ${\widetilde{\mathbf{u}}_{\E_0}}$ can be found as ${\widetilde{\mathbf{u}}_{\E_0}}= \lim_{n\rightarrow \infty} {\mathbf{u}}_{\E_0}^{(n)}$ using the following iterative procedure representing a discretization of a gradient flow projected on $\mathcal{S}_{\E_0}$ $$\begin{aligned}
{\mathbf{u}}_{\E_0}^{(n+1)} & =  \mathbb{P}_{\mathcal{S}_{\E_0}}\left(\;{\mathbf{u}}^{(n)}_{\E_0} + \tau_n \nabla\R\left({\mathbf{u}}^{(n)}_{\E_0}\right)\;\right), \\ 
{\mathbf{u}}_{\E_0}^{(1)} & =  {\mathbf{u}}^0,
\end{aligned}
\label{eq:desc}$$ where ${\mathbf{u}}^{(n)}_{\E_0}$ is an approximation of the maximizer obtained at the $n$-th iteration, ${\mathbf{u}}^0$ is the initial guess and $\tau_n$ is the length of the step in the direction of the gradient. It is ensured that the maximizers ${\widetilde{\mathbf{u}}_{\E_0}}$ obtained for different values of $\E_0$ lie on the same [maximizing]{} branch by using the continuation approach, where the maximizer ${\widetilde{\mathbf{u}}_{\E_0}}$ is [employed]{} as the initial guess ${\mathbf{u}}^0$ to compute $\widetilde{\mathbf{u}}_{\E_0+\Delta\E}$ at the next enstrophy level for some sufficiently small $\Delta\E > 0$. As will be demonstrated in §\[sec:3D\_InstOpt\_E0to0\], in the limit $\E_0 \rightarrow 0$ optimization problem \[pb:maxdEdt\_E\] admits a discrete family of closed-form solutions and each of these vortex states is the limiting (initial) member $\widetilde{\mathbf{u}}_{0}$ of the corresponding [maximizing]{} branch. As such, these limiting extreme vortex states are used as the initial guesses ${\mathbf{u}}^0$ for the calculation of $\widetilde{\mathbf{u}}_{\Delta\E}$, i.e., they serve as “seeds” for the calculation of an entire [maximizing]{} branch (as discussed in §\[sec:discuss\], while there exist alternatives to the continuation approach, this technique [in fact results]{} in the fastest convergence of iterations and also ensures that all computed extreme vortex states lie on a single branch). The procedure outlined above is summarized as Algorithm \[alg:optimAlg\], [whereas all]{} details are presented below.

set $\E_0 = 0$ set ${\widetilde{\mathbf{u}}_{\E_0}}= \widetilde{\mathbf{u}}_{0}$ ${\mathbf{u}}_{\E_0}^{(0)} = {\widetilde{\mathbf{u}}_{\E_0}}$ $\E_0 = \E_0 + \Delta \E$ $n = 0$ compute $\R_0 = \R\left({\mathbf{u}}_{\E_0}^{(0)}\right)$

compute the $L_2$ gradient $\nabla^{L_2}\R\left({\mathbf{u}}_{\E_0}^{(n)}\right)$, see equation

compute the Sobolev gradient $\nabla\R\left({\mathbf{u}}_{\E_0}^{(n)}\right)$, see equation

compute the step size $\tau_n$, see equation

set ${\mathbf{u}}_{\E_0}^{(n+1)} = \mathbb{P}_{\mathcal{S}_{\E_0}}\left(\;{\mathbf{u}}_{\E_0}^{(n)} + \tau_n \nabla\R\left({\mathbf{u}}_{\E_0}^{(n)}\right)\;\right)$, see equations –

set $\R_1 = \R\left({\mathbf{u}}_{\E_0}^{(n+1)}\right)$

compute the `relative error` $ = (\R_1 - \R_0)/\R_0$

set $\R_0 = \R_1$

set $n=n+1$

A key step of Algorithm \[alg:optimAlg\] is the evaluation of the gradient $\nabla\R({\mathbf{u}})$ of the objective functional $\R({\mathbf{u}})$, cf. , representing its (infinite-dimensional) sensitivity to perturbations of the velocity field ${\mathbf{u}}$, and it is essential that the gradient be characterized by the required regularity, namely, $\nabla\R({\mathbf{u}}) \in H^2(\Omega)$. This is, in fact, guaranteed by the Riesz representation theorem [@l69] applicable because the Gâteaux differential $\R'({\mathbf{u}};\cdot) : H_0^2(\Omega) \rightarrow \RR$, defined as $\R'({\mathbf{u}};{\mathbf{u}}') := \lim_{\epsilon \rightarrow 0}
\epsilon^{-1}\left[\R({\mathbf{u}}+\epsilon {\mathbf{u}}') - \R({\mathbf{u}})\right]$ for some perturbation ${\mathbf{u}}' \in H_0^2(\Omega)$, is a bounded linear functional on $H_0^2(\Omega)$. The Gâteaux differential can be computed directly to give $$\R'({\mathbf{u}};{\mathbf{u}}') = \int_{\Omega}\left[{\mathbf{u}}'\cdot\bnabla{\mathbf{u}}\cdot{\Delta}{\mathbf{u}}+ 
{\mathbf{u}}\cdot\bnabla{\mathbf{u}}'\cdot{\Delta}{\mathbf{u}}+ 
{\mathbf{u}}\cdot\bnabla{\mathbf{u}}\cdot{\Delta}{\mathbf{u}}' \right]\,d{\mathbf{x}}-2\nu\int_{\Omega}{\Delta}^2{\mathbf{u}}\cdot{\mathbf{u}}'\,d{\mathbf{x}}\label{eq:dR}$$ from which, by the Riesz representation theorem, we obtain $$\R'({\mathbf{u}};{\mathbf{u}}') 
= \Big\langle \nabla\R({\mathbf{u}}), {\mathbf{u}}' \Big\rangle_{H^2(\Omega)}
= \Big\langle \nabla^{L_2}\R({\mathbf{u}}), {\mathbf{u}}' \Big\rangle_{L_2(\Omega)}
\label{eq:riesz}$$ with the Riesz representers $\nabla\R({\mathbf{u}})$ and $\nabla^{L_2}\R({\mathbf{u}})$ being the gradients computed with respect to the $H^2$ and $L_2$ topology, respectively, and the inner products defined in . We remark that, while the $H^2$ gradient is used exclusively in the actual computations, cf. , the $L_2$ gradient is computed first as an intermediate step. Identifying the Gâteaux differential with the $L_2$ inner product and performing integration by parts yields $$\nabla^{L_2}\R({\mathbf{u}}) = {\Delta}\left( {\mathbf{u}}\cdot\bnabla{\mathbf{u}}\right) + (\bnabla{\mathbf{u}})^T{\Delta}{\mathbf{u}}- 
{\mathbf{u}}\cdot\bnabla({\Delta}{\mathbf{u}}) - 2\nu{\Delta}^2{\mathbf{u}}.
\label{eq:gradRL2}$$ Similarly, identifying the Gâteaux differential with the $H^2$ inner product , integrating by parts and using , we obtain the required $H^2$ gradient $\nabla\R$ as a solution of the following elliptic boundary-value problem $$\begin{aligned}
&\left[ {\operatorname{Id}}\, - \,\ell_1^2 \,\Delta + \,\ell_2^4 \,\Delta^2 \right] \nabla\R
= \nabla^{L_2} \R  \qquad \text{in} \ \Omega, \\
& \text{Periodic Boundary Conditions}.
\end{aligned}
\label{eq:gradRH2}$$ The gradient fields $\nabla^{L_2}\R({\mathbf{u}})$ and $\nabla\R({\mathbf{u}})$ can be interpreted as infinite-dimensional sensitivities of the objective function $\R({\mathbf{u}})$, cf. , with respect to perturbations of the field ${\mathbf{u}}$. While these two gradients may point towards the same local maximizer, they represent distinct “directions”, since they are defined with respect to different topologies ($L_2$ vs. $H^2$). As shown by @pbh04, extraction of gradients in spaces of smoother functions such as $H^2(\Omega)$ can be interpreted as low-pass filtering of the $L_2$ gradients with parameters $\ell_1$ and $\ell_2$ acting as cut-off length-scales and the choice of their numerical values will be discussed in §\[sec:3D\_InstOpt\_E\].

The step size $\tau_n$ in algorithm is computed as $$\label{eq:tau_n}
\tau_n = \mathop{{\operatorname{argmax}}}_{\tau>0} \left\{ \R\left[\mathbb{P}_{\mathcal{S}_{\E_0}}
\left( \;{\mathbf{u}}^{(n)} + \tau\,\nabla\R({\mathbf{u}}^{(n)}) \;\right)\right] \right\}$$ which is done using a suitable derivative-free line-search algorithm [@r06]. Equation can be interpreted as a modification of a standard line search method where the optimization is performed following an arc (a geodesic) lying on the constraint manifold $\mathcal{S}_{\E_0}$, rather than a straight line. This approach was already successfully employed to solve similar problems in @ap11a [@ap13a].

It ought to be emphasized here that the approach presented above in which the projections – and gradients – are obtained based on the infinite-dimensional (continuous) formulation to be discretized only at the final stage is fundamentally different from the method employed in the original study by [@ld08] in which the optimization problem was solved in a fully discrete setting (the two approaches are referred to as “optimize-then-discretize” and “discretize-then-optimize”, respectively, cf. [@g03]). A practical advantage of the continuous (“optimize-then-discretize”) formulation used in the present work is that the expressions representing the sensitivity of the objective functional $\R$, i.e. the gradients $\nabla^{L_2}\R$ and $\nabla\R$, are independent of the specific discretization approach chosen to evaluate them. This should be contrasted with the discrete (“discretize-then-optimize”) formulation, where a change of the discretization method would require rederivation of the gradient expressions. In addition, the continuous formulation allows us to strictly enforce the regularity of maximizers required in problem \[pb:maxdEdt\_E\]. Finally and perhaps most importantly, the continuous formulation of the maximization problem makes it possible to obtain elegant closed-form solutions of the problem in the limit $\E_0 \rightarrow 0$, which is done in § \[sec:3D\_InstOpt\_E0to0\] below. These analytical solutions will then be used in §\[sec:3D\_InstOpt\_E\] to guide the computation of maximizing branches by numerically solving problem \[pb:maxdEdt\_E\] for a broad range of $\E_0$, as outlined in Algorithm \[alg:optimAlg\].

Extreme Vortex States in the Limit $\E_0 \to 0$ {#sec:3D_InstOpt_E0to0}
===============================================

It is possible to find analytic solutions to problem \[pb:maxdEdt\_E\] in the limit $\E_0 \to 0$ using perturbation methods. To simplify the notation, in this section we will drop the subscript $\E_0$ when referring to the optimal field. The Euler-Lagrange system representing the first-order optimality conditions in optimization problem \[pb:maxdEdt\_E\] is given by [@l69]

\[eq:KKT\_E\] $$\begin{aligned}
{\mathcal{B}}({\widetilde{\mathbf{u}}},{\widetilde{\mathbf{u}}}) - 2\nu{\Delta}^2{\widetilde{\mathbf{u}}}- \lambda{\Delta}{\widetilde{\mathbf{u}}}- \bnabla q & = 0 \qquad\mbox{in}\,\,\Omega , \label{eq:KKT_E_gradR}\\
\nabla\cdot{\widetilde{\mathbf{u}}}& = 0 \qquad\mbox{in}\,\,\Omega , \label{eq:KKT_E_divConstr}\\
\E({\widetilde{\mathbf{u}}}) - \E_0 & = 0, \label{eq:KKT_E_E0Constr}\end{aligned}$$

where $\lambda\in\mathbb{R}$ and $q:\Omega\to\mathbb{R}$ are the Lagrange multipliers associated with the constraints defining the manifold ${\mathcal{S}_{\E_0}}$, and ${\mathcal{B}}({\mathbf{u}},{\mathbf{v}})$, given by $${\mathcal{B}}({\mathbf{u}},{\mathbf{v}}) :=  {\Delta}\left( {\mathbf{u}}\cdot\bnabla{\mathbf{v}}\right) + (\bnabla{\mathbf{u}})^T{\Delta}{\mathbf{v}}- 
{\mathbf{u}}\cdot\bnabla({\Delta}{\mathbf{v}}),$$ is the bilinear form from equation . Using the formal series expansions with $\alpha > 0$

\[eq:series3D\] $$\begin{aligned}
{\widetilde{\mathbf{u}}}& = {\mathbf{u}}_0 + \E_0^{\alpha}{\mathbf{u}}_1 + \E_0^{2\alpha}{\mathbf{u}}_2 + \ldots, \\
\lambda & = \lambda_0 + \E_0^{\alpha}\lambda_1 + \E_0^{2\alpha}\lambda_2 + \ldots, \\
q & = q_0 + \E_0^{\alpha}q_1 + \E_0^{2\alpha}q_2 + \ldots\end{aligned}$$

in and collecting terms proportional to different powers of $\E_0^{\alpha}$, it follows from that, at every order $m=1,2,\dots$ in $\E_0^{\alpha}$, we have $$\E_0^{m\alpha}: \qquad\sum_{j=0}^m {\mathcal{B}}({\mathbf{u}}_j,{\mathbf{u}}_{m-j}) - 2\nu{\Delta}^2{\mathbf{u}}_m -
\sum_{j=0}^m\lambda_j{\Delta}{\mathbf{u}}_{m-j} - \nabla q_m = 0 \quad\mbox{in}\,\,\Omega.$$ Similarly, equation leads to $$\label{eq:Incompressible_Uk} 
\nabla\cdot{\mathbf{u}}_m = 0 \quad\mbox{in}\,\,\Omega$$ at every order $m$ in $\E_0^{\alpha}$. It then follows from equation that $$\begin{aligned}
\E({\mathbf{u}}) & = & \E({\mathbf{u}}_0) -\big\langle{\mathbf{u}}_0,{\Delta}{\mathbf{u}}_1\big\rangle_{L_2}\E_0^{\alpha} + 
                    \left[ \E({\mathbf{u}}_1) - \big\langle{\mathbf{u}}_0,{\Delta}{\mathbf{u}}_2\big\rangle_{L_2}\right]\E_0^{2\alpha} + \ldots \\
              & = & \E_0,\end{aligned}$$ which, for $\alpha \neq 0$, forces $\E({\mathbf{u}}_0) = 0$. Hence, ${\mathbf{u}}_0
\equiv 0$, $\alpha = 1/2$ and $\E({\mathbf{u}}_1) = 1$. The systems at orders $\E_0^{1/2}$ and $\E_0^1$ are given by:

\[eq:maxdEdt\_Asympt\_1\] $$\begin{aligned}
\E^{1/2}_0:\quad\qquad\qquad\qquad\qquad\qquad 2\nu{\Delta}^2{\mathbf{u}}_1 + \lambda_0{\Delta}{\mathbf{u}}_1 + 
\nabla q_1 & = 0 \quad\mbox{in}\,\,\Omega , \label{eq:maxdEdt_Asympt_1_PDE}\\
\nabla\cdot{\mathbf{u}}_1 & = 0 \quad\mbox{in}\,\,\Omega , \label{eq:maxdEdt_Asympt_1_Div0Constr} \\
\E({\mathbf{u}}_1) & = 1, \label{eq:maxdEdt_Asympt_1_E0Constr} \end{aligned}$$

\[eq:maxdEdt\_Asympt\_2\] $$\begin{aligned}
\E_0:\qquad 2\nu{\Delta}^2{\mathbf{u}}_2 + \lambda_0{\Delta}{\mathbf{u}}_2  + \nabla q_2 - {\mathcal{B}}({\mathbf{u}}_1,{\mathbf{u}}_1) + 
\lambda_1{\Delta}{\mathbf{u}}_1 & = 0 \quad\mbox{in}\,\,\Omega, \\
\nabla\cdot{\mathbf{u}}_2 & = 0 \quad\mbox{in}\,\,\Omega , \\
\langle {\Delta}{\mathbf{u}}_1, {\mathbf{u}}_2 \rangle_{L_2} & = 0, \end{aligned}$$

where the fact that ${\mathcal{B}}({\mathbf{u}}_0,{\mathbf{u}}_j) = 0$ for all $j$ has been used. While continuing this process to larger values of $m$ may lead to some interesting insights, for the purpose of this investigation it is sufficient to truncate expansions at the order $\O(\E_0)$. The corresponding approximation of the objective functional then becomes $$\label{eq:R03D}
\R({\widetilde{\mathbf{u}}}) = - \nu\E_0\int_{\Omega} \left| {\Delta}{\mathbf{u}}_1 \right|^2 \, d{\mathbf{x}}+ \O(\E_0^{3/2}).$$ It is worth noting that, in the light of relation , the maximum rate of growth of enstrophy in the limit of small $\E_0$ is in fact negative, meaning that, for sufficiently small $\E_0$, the enstrophy itself is a decreasing function for all times. This observation is consistent with the small-data regularity result discussed in Introduction.

As regards problem defining the triplet $\{{\mathbf{u}}_1, q_1,\lambda_0\}$, taking the divergence of equation and using the condition $\nabla\cdot{\mathbf{u}}_1=0$ leads to the Laplace equation ${\Delta}q_1 =
0$ in $\Omega$. Since for zero-mean functions defined on $\Omega$, ${\operatorname{Ker}}({\Delta}) = \{ 0 \}$, it follows that $q_1
\equiv 0$ and equation is reduced to the eigenvalue problem $$\label{eq:maxdEdt_smallE0_eig}
2\nu{\Delta}{\mathbf{u}}_1 + \lambda_0{\mathbf{u}}_1 = 0,$$ with ${\mathbf{u}}_1$ satisfying the incompressibility condition . Direct calculation using equation and condition leads to an asymptotic expression for the objective functional in the limit of small enstrophy $$\label{eq:R0_lambda0_3D}
\R({\widetilde{\mathbf{u}}}) \approx - \lambda_0\E_0.$$ Solutions to the eigenvalue problem in equation can be found using the Fourier expansion of ${\mathbf{u}}_1$ given as (with hats denoting Fourier coefficients) $${\mathbf{u}}_1({\mathbf{x}}) = \sum_{{\mathbf{k}}\in{\mathcal{W}}}\widehat{{\mathbf{u}}}_1({\mathbf{k}}){\textrm{e}^{2\pi i{\mathbf{k}}\cdot{\mathbf{x}}}},$$ where ${\mathcal{W}}\subseteq\mathbb{Z}^3$ is a set of wavevectors ${\mathbf{k}}$ for which $\widehat{{\mathbf{u}}}_1({\mathbf{k}}) \neq 0$. The eigenvalue problem then becomes $$\begin{aligned}
\left[-2\nu(2\pi)^2|{\mathbf{k}}|^2 + \lambda_0\right]\widehat{{\mathbf{u}}}_1({\mathbf{k}}) & = & 0 \qquad\forall\,{\mathbf{k}}\in{\mathcal{W}}, \\
\widehat{{\mathbf{u}}}_1({\mathbf{k}})\cdot{\mathbf{k}}& = & 0 \qquad\forall\,{\mathbf{k}}\in{\mathcal{W}},\end{aligned}$$ with solutions obtained by choosing, for any $k\in\mathbb{Z}\setminus\{0\}$, a set of wavevectors with the following structure $$\label{eq:defW}
{\mathcal{W}}_k = \left\{{\mathbf{k}}\in\mathbb{Z}^3\colon|{\mathbf{k}}|^2 = k \right\}$$ and $\widehat{{\mathbf{u}}}_1({\mathbf{k}})$ with an appropriate form satisfying the incompressibility condition $\widehat{{\mathbf{u}}}_1\cdot{\mathbf{k}}= 0$. For the solutions to equation constructed in such manner it then follows that $\lambda_0 = 2\nu(2\pi)^2|{\mathbf{k}}|^2$ and the optimal asymptotic value of $\R$ is given by $$\label{eq:R0_kvec_3D}
\R({\widetilde{\mathbf{u}}}) \approx - 8\pi^2\nu|{\mathbf{k}}|^2\E_0.$$

Since the fields ${\mathbf{u}}_1$ are real-valued, their Fourier modes must satisfy $\widehat{{\mathbf{u}}}_1(-{\mathbf{k}}) =
\overline{\widehat{{\mathbf{u}}}_1({\mathbf{k}})}$, where $\overline{z}$ denotes the complex conjugate (C.C.) of $z\in\mathbb{C}$. Depending on the choice of ${\mathcal{W}}_k$, a number of different solutions of can be constructed and below we focus on the following three most relevant cases characterized by the largest values of $\R({\widetilde{\mathbf{u}}})$:

1.  ${\mathcal{W}}_1 = \{ {\mathbf{k}}_1, {\mathbf{k}}_2, {\mathbf{k}}_3, -{\mathbf{k}}_1, -{\mathbf{k}}_2,
      -{\mathbf{k}}_3 \}$, where ${\mathbf{k}}_i = \mathbf{e}_i$, $i=1,2,3$, is the $i^{\textrm{th}}$ unit vector of the canonical basis of $\mathbb{R}^3$; the most general solution can then be constructed as $$\label{eq:uvec_3D_k1}
    {\mathbf{u}}_1({\mathbf{x}}) = \mathbf{A}{\textrm{e}^{2\pi i{\mathbf{k}}_1\cdot{\mathbf{x}}}} + 
                     \mathbf{B}{\textrm{e}^{2\pi i{\mathbf{k}}_2\cdot{\mathbf{x}}}} +
                     \mathbf{C}{\textrm{e}^{2\pi i{\mathbf{k}}_3\cdot{\mathbf{x}}}} + \textrm{C.C.}$$ with the complex-valued constant vectors $\mathbf{A} = [0,A_2,A_3]$, $\mathbf{B} = [B_1,0,B_3]$ and $\mathbf{C} = [C_1,C_2,0]$ suitably chosen so that $\E({\mathbf{u}}_1) = 1$; hereafter we will use the values $A_2 =
    A_3 = \ldots = C_2 = 1/(48\pi^2)$; it follows that $|{\mathbf{k}}|^2 =
    1$ $\forall\,\,{\mathbf{k}}\in{\mathcal{W}}_1$, and the optimal asymptotic value of $\R$ obtained from equation is given by $$\label{eq:R0_kvec_3D_k1}
      \R({\widetilde{\mathbf{u}}}) \approx - 8\pi^2\nu\E_0,$$ \[c1\]

2.  ${\mathcal{W}}_2 = {\mathcal{W}}\cup (-{\mathcal{W}})$, where $-{\mathcal{W}}$ denotes the set whose elements are the negatives of the elements of set ${\mathcal{W}}$, for ${\mathcal{W}}= \{
      {\mathbf{k}}_1 + {\mathbf{k}}_2, {\mathbf{k}}_1 - {\mathbf{k}}_2, {\mathbf{k}}_1 + {\mathbf{k}}_3, {\mathbf{k}}_1 -
      {\mathbf{k}}_3, {\mathbf{k}}_2 + {\mathbf{k}}_3, {\mathbf{k}}_2 - {\mathbf{k}}_3 \}$; the most general solution can be then constructed as $$\begin{aligned}
    {\mathbf{u}}_1({\mathbf{x}}) & = & \mathbf{A}{\textrm{e}^{2\pi i[1,1,0]\cdot{\mathbf{x}}}} + 
                         \mathbf{B}{\textrm{e}^{2\pi i[1,-1,0]\cdot{\mathbf{x}}}} +
                         \mathbf{C}{\textrm{e}^{2\pi i[1,0,1]\cdot{\mathbf{x}}}} + \nonumber \\
                   &   & \mathbf{D}{\textrm{e}^{2\pi i[1,0,-1]\cdot{\mathbf{x}}}} +
                         \mathbf{E}{\textrm{e}^{2\pi i[0,1,1]\cdot{\mathbf{x}}}} + 
                         \mathbf{F}{\textrm{e}^{2\pi i[0,1,-1]\cdot{\mathbf{x}}}} + \textrm{C.C.} 
    \label{eq:uvec_3D_k2}\end{aligned}$$ with the constants $\mathbf{A},\mathbf{B},\ldots,\mathbf{F}\in\mathbb{C}^3$ suitably chosen so that $\mathbf{A}\cdot[1,1,0] = 0$, $\mathbf{B}\cdot[1,-1,0]
    = 0 ,\ldots,\mathbf{F}\cdot[0,1,-1] = 0$, which ensures that incompressibility condition is satisfied, and that $\E({\mathbf{u}}_1) = 1$; in this case, $|{\mathbf{k}}|^2 =
    2$, $\forall\,{\mathbf{k}}\in{\mathcal{W}}_2$, and the optimal asymptotic value of $\R$ is $$\label{eq:R0_kvec_3D_k2}
    \R({\widetilde{\mathbf{u}}}) \approx - 16\pi^2\nu\E_0,$$ \[c2\]

3.  ${\mathcal{W}}_3 = {\mathcal{W}}\cup (-{\mathcal{W}})$ for ${\mathcal{W}}= \{
      {\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3,-{\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3,{\mathbf{k}}_1-{\mathbf{k}}_2+{\mathbf{k}}_3,{\mathbf{k}}_1+{\mathbf{k}}_2-{\mathbf{k}}_3
      \}$; the most general solution can then be constructed as $$\begin{aligned}
    {\mathbf{u}}_1({\mathbf{x}}) & = & \mathbf{A}{\textrm{e}^{2\pi i[1,1,1]\cdot{\mathbf{x}}}} + 
                         \mathbf{B}{\textrm{e}^{2\pi i[-1,1,1]\cdot{\mathbf{x}}}} + \nonumber \\
                   &   & \mathbf{C}{\textrm{e}^{2\pi i[1,-1,1]\cdot{\mathbf{x}}}} + 
                         \mathbf{D}{\textrm{e}^{2\pi i[1,1,-1]\cdot{\mathbf{x}}}} + \textrm{C.C.} 
    \label{eq:uvec_3D_k3}\end{aligned}$$ with the constants $\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}\in\mathbb{C}^3$ suitably chosen so that $\mathbf{A}\cdot[1,1,1] = 0$, $\mathbf{B}\cdot[-1,1,1]
    = 0$, $\mathbf{C}\cdot[1,-1,1] = 0$ and $\mathbf{D}\cdot[1,1,-1] = 0$, which ensures that incompressibility condition is satisfied, and that $\E({\mathbf{u}}_1) = 1$; in this case, $|{\mathbf{k}}|^2 =3$, $\forall\,{\mathbf{k}}\in{\mathcal{W}}_3$, and the optimal asymptotic value of $\R$ is $$\label{eq:R0_kvec_3D_k3}
    \R({\widetilde{\mathbf{u}}}) \approx - 24\pi^2\nu\E_0.$$ \[c3\]

The three constructions of the extremal field ${\mathbf{u}}_1$ given in , and are all defined up to arbitrary shifts in all three directions, reflections with respect to different planes and rotations by angles which are multiples of $\pi / 2$ about the different axes. As a result of this nonuniqueness, there is some freedom in choosing the constants $\mathbf{A},\ldots,\mathbf{F}$. Given that the optimal asymptotic value of $\R$ depends exclusively on the wavevector magnitude $|{\mathbf{k}}|$, cf. , any combination of constants $\mathbf{A},\ldots,\mathbf{F}$ will produce the [*same*]{} optimal rate of growth of enstrophy. Thus, to fix attention, in our analysis we will set $\mathbf{A}=\mathbf{B}=\mathbf{C}$ in case (i), $\mathbf{A}=\mathbf{B}=\ldots=\mathbf{F}$ in case (\[c2\]) and $\mathbf{A}=\ldots=\mathbf{D}$ in case (\[c3\]). With these choices, the contribution from each component of the field ${\mathbf{u}}_1$ to the total enstrophy is the same. The maximum (i.e., least negative) value of $\R$ can be thus obtained by choosing the smallest possible $|{\mathbf{k}}|^2$. This maximum is achieved in case (\[c1\]) with the wavevectors ${\mathbf{k}}_1 = [1,0,0]$, ${\mathbf{k}}_2 = [0,1,0]$, ${\mathbf{k}}_3 =
[0,0,1]$, and $-{\mathbf{k}}_1$, $-{\mathbf{k}}_2$ and $-{\mathbf{k}}_3$, for which $|{\mathbf{k}}|^2 = 1$. Because of this maximization property, this is the field we will focus on in our analysis in §\[sec:3D\_InstOpt\_E\] and §\[sec:timeEvolution\].

The three fields constructed in , and are visualized in figure \[fig:maxdEdt\_vortexCells\]. This analysis is performed using the level sets $\Gamma_{s}(F)\subset\Omega$ defined as $$\label{eq:levelSets}
\Gamma_{s}(F) :=  \{{\mathbf{x}}\in\Omega : F({\mathbf{x}}) = s \},$$ for a suitable function $F:\Omega\to\mathbb{R}$. In figures \[fig:maxdEdt\_vortexCells\](a–c) we choose $F({\mathbf{x}}) =
  |{\bnabla\times}{\mathbf{u}}_1|({\mathbf{x}})$ with $s = 0.95||{\bnabla\times}{\mathbf{u}}_1||_{L_\infty}$. To complement this information, in figures \[fig:maxdEdt\_vortexCells\](d–f) we also plot the isosurfaces and cross-sectional distributions of the $x_1$ component of the field ${\mathbf{u}}_1$.

The fields shown in figure \[fig:maxdEdt\_vortexCells\] reveal interesting patterns involving well-defined “vortex cells”. More specifically, we see that in case (\[c1\]), given by equation and shown in figures \[fig:maxdEdt\_vortexCells\](a,d), the vortex cells are staggered with respect to the orientation of the cubic domain $\Omega$ in all three planes, whereas in case (\[c3\]), given by equation and shown in figures \[fig:maxdEdt\_vortexCells\](c,f), the vortex cells are aligned with the domain $\Omega$ in all three planes. On the other hand, in case (\[c2\]), given by equation and shown in figures \[fig:maxdEdt\_vortexCells\](b,e), the vortex cells are staggered in one plane and aligned in another with the arrangement in the third plane resulting from the arrangement in the first two. These geometric properties are also reflected in the $x_1$-component of the field ${\mathbf{u}}_1$ which is independent of $x_1$ in cases (\[c1\]) and (\[c2\]), but exhibits, respectively, a staggered and aligned arrangement of the cells in the $y-z$ plane in these two cases. In case (\[c3\]) the cells exhibit an aligned arrangement in all three planes. The geometric properties of the extreme vortex states obtained in the limit $\E_0 \rightarrow 0$ are summarized in Table \[tab:E0\]. We remark that an analogous structure of the optimal fields, featuring aligned and staggered arrangements of vortex cells in the limiting case, was also discovered by [@ap13a] in their study of the maximum palinstrophy growth in 2D. While due to a smaller spatial dimension only two optimal solutions were found in that study, the one characterized by the staggered arrangement also lead to a larger (less negative) rate of palinstrophy production.

[l|c|c|c|c]{}

[Case]{}\

&

[Formula for the velocity field]{}

&

[Arrangement of cells in $y-z$ plane]{}

&

[Dependence of $x_1$ component of ${\mathbf{u}}_1$ on $x_1$]{}

&

[Remarks]{}

\

(\[c1\])\

&

&

staggered

&

uniform

&

staggered ABC flow

\

(\[c2\])\

&

&

aligned

&

uniform

&

aligned ABC flow

\

(\[c3\])\

&

&

aligned

&

cell-like

&

Taylor-Green vortex

\

It is also worth mentioning that the initial data for two well-known flows, namely, the Arnold-Beltrami-Childress (ABC) flow [@mb02] and the Taylor-Green flow [@tg37], are in fact particular instances of the optimal field ${\mathbf{u}}_1$ corresponding to, respectively, cases (\[c1\]) and (\[c3\]). Following the notation of @df86, general ABC flows are characterized by the following velocity field $$\label{eq:ABC_flow}
\begin{array}{r@{\,\,}c@{\,\,}l}
u_1(x_1,x_2,x_3) & = & A'\sin(2\pi x_3) + C'\cos(2\pi x_2), \\
u_2(x_1,x_2,x_3) & = & B'\sin(2\pi x_1) + A'\cos(2\pi x_3), \\
u_3(x_1,x_2,x_3) & = & C'\sin(2\pi x_2) + B'\cos(2\pi x_1),
\end{array}$$ where $A'$, $B'$ and $C'$ real constants. The vector field in equation can be obtained from equation by choosing $\mathbf{A} = (B'/2)[0,-i,1]$, $\mathbf{B} = (C'/2)[1,0,-i]$ and $\mathbf{C} = (A'/2)[-i, 1,0]$. By analogy, we will refer to the state described by as the “aligned ABC flow”. Likewise, the well-known Taylor-Green vortex can be obtained as a particular instance of the field ${\mathbf{u}}_1$ from equation again using a suitable choice of the constants $\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}$. Traditionally, the velocity field ${\mathbf{u}}= [u_1,u_2,u_3]$ characterizing the Taylor-Green vortex is defined as [@bmonmu83] $$\label{eq:TaylorGreen_vortex}
\begin{array}{r@{\,\,}c@{\,\,}l}
u_1(x_1,x_2,x_3) & = & \gamma_1\sin(2\pi x_1)\cos(2\pi x_2)\cos(2\pi x_3), \\
u_2(x_1,x_2,x_3) & = & \gamma_2\cos(2\pi x_1)\sin(2\pi x_2)\cos(2\pi x_3), \\
u_3(x_1,x_2,x_3) & = & \gamma_3\cos(2\pi x_1)\cos(2\pi x_2)\sin(2\pi x_3), \\
 0 & = & \gamma_1+\gamma_2+\gamma_3,  
\end{array}$$ for $\gamma_1,\gamma_2,\gamma_3\in\mathbb{R}$. For given values of $\gamma_1$,$ \gamma_2$ and $\gamma_3$ in , the corresponding constants $\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}$ in can be found by separating them into their real and imaginary parts denoted, respectively, ${\mathbf{A}_{\textrm{Re}}},{\mathbf{B}_{\textrm{Re}}},{\mathbf{C}_{\textrm{Re}}},{\mathbf{D}_{\textrm{Re}}}$ and ${\mathbf{A}_{\textrm{Im}}},{\mathbf{B}_{\textrm{Im}}},{\mathbf{C}_{\textrm{Im}}},{\mathbf{D}_{\textrm{Im}}}$. Then, after choosing $${\mathbf{A}_{\textrm{Re}}}={\mathbf{B}_{\textrm{Re}}}={\mathbf{C}_{\textrm{Re}}}={\mathbf{D}_{\textrm{Re}}}=\mathbf{0} = [0,0,0],$$ the imaginary parts can be determined by solving the following system of linear equations $$2\left[
\begin{array}{cccc}
I_3 & -I_3 & -I_3 & -I_3 \\ 
-I_3 & I_3 & -I_3 & -I_3 \\ 
-I_3 & -I_3 & I_3 & -I_3 \\ 
-I_3 & -I_3 & -I_3 & I_3   
\end{array} 
\right]
\left[
\begin{array}{c}
{\mathbf{A}_{\textrm{Im}}} \\
{\mathbf{B}_{\textrm{Im}}} \\
{\mathbf{C}_{\textrm{Im}}} \\
{\mathbf{D}_{\textrm{Im}}}
\end{array}
\right] = 
\left[
\begin{array}{c}
\mathbf{0} \\
\gamma_1\mathbf{e}_1 \\
\gamma_2\mathbf{e}_2 \\
\gamma_3\mathbf{e}_3
\end{array}
\right],$$ where $I_3$ is the $3 \times 3$ identity matrix. The values of ${\mathbf{A}_{\textrm{Im}}},\ldots,{\mathbf{D}_{\textrm{Im}}}$ are thus given by $${\mathbf{A}_{\textrm{Im}}} = -\frac{1}{8}\left[ \begin{array}{c}\gamma_1 \\ \gamma_2 \\ \gamma_3 \end{array} \right],\, 
{\mathbf{B}_{\textrm{Im}}} = -\frac{1}{8}\left[ \begin{array}{c}-\gamma_1 \\ \gamma_2 \\ \gamma_3 \end{array} \right],\,
{\mathbf{C}_{\textrm{Im}}} = -\frac{1}{8}\left[ \begin{array}{c}\gamma_1 \\ -\gamma_2 \\ \gamma_3 \end{array} \right],\,
{\mathbf{D}_{\textrm{Im}}} = -\frac{1}{8}\left[ \begin{array}{c}\gamma_1 \\ \gamma_2 \\ -\gamma_3 \end{array} \right].$$ A typical choice of the parameters used in the numerical studies performed by [@bmonmu83] and [@b91] is $\gamma_1 = - \gamma_2
= 1$ and $\gamma_3 = 0$.

We remark that the Taylor-Green vortex has been employed as the initial data in a number of studies aimed at triggering singular behaviour in both the Euler and Navier-Stokes systems [@tg37; @bmonmu83; @b91; @bb12]. It is therefore interesting to note that it arises in the limit $\E_0 \rightarrow 0$ as one of the extreme vortex states in the variational formulation considered in the present study. It should be emphasized, however, that out of the three optimal states identified above (see Table \[tab:E0\]), the Taylor-Green vortex is characterized by the smallest (i.e., the most negative) instantaneous rate of enstrophy production $d\E/dt$. On the other hand, we are not aware of any prior studies involving ABC flows in the context of extreme behaviour and potential singularity formation. The time evolution corresponding to these states and some other initial data will be analyzed in detail in §\[sec:timeEvolution\].

Extreme Vortex States with Finite $\E_0$ {#sec:3D_InstOpt_E}
========================================

In this section we analyze the optimal vortex states ${\widetilde{\mathbf{u}}_{\E_0}}$ obtained for finite values of the enstrophy in which we extend the results obtained in the seminal study by [@ld08]. As was also the case in the analogous study in 2D [@ap13a], there is a distinct branch of extreme states ${\widetilde{\mathbf{u}}_{\E_0}}$ parameterized by the enstrophy $\E_0$ and corresponding to each of the three limiting states discussed in §\[sec:3D\_InstOpt\_E0to0\] (cf. figure \[fig:maxdEdt\_vortexCells\] and Table \[tab:E0\]). Each of these branches is computed using the continuation approach outlined in Algorithm \[alg:optimAlg\]. As a key element of the gradient-based maximization technique , the gradient expressions – are approximated pseudo-spectrally using standard dealiasing of the nonlinear terms and with resolutions varying from $128^3$ in the low-enstrophy cases to $512^3$ in the high-enstrophy cases, which necessitated a massively parallel implementation using the Message Passing Interface (MPI). As regards the computation of the Sobolev $H^2$ gradients, cf. , we set $\ell_1 = 0$, whereas the second parameter $\ell_2$ was adjusted during the optimization iterations and was chosen so that $\ell_2 \in [\ell_{\min},
\ell_{\max}]$, where $\ell_{\min}$ is the length scale associated with the spatial resolution $N$ used for computations and $\ell_{\max}$ is the characteristic length scale of the domain $\Omega$, that is, $\ell_{\min} \sim \O( 1/N) $ and $\ell_{\max} \sim \O(1)$. We remark that, given the equivalence of the inner products corresponding to different values of $\ell_1$ and $\ell_2$ (as long as $\ell_2 \neq 0$), these choices do not affect the maximizers found, but only how rapidly they are approached by iterations . For further details concerning the computational approach we refer the reader to the dissertation by [@a14]. As was the case in the analogous 2D problem studied by @ap13a, the largest instantaneous growth of enstrophy is produced by the states with vortex cells staggered in all planes, cf. case (\[c1\]) in Table \[tab:E0\]. Therefore, in our analysis we will focus exclusively on this branch of extreme vortex states which has been computed for $\E_0 \in [10^{-3},2\times10^2]$.

The optimal instantaneous rate of growth of enstrophy $\R_{\E_0} =
\R({\widetilde{\mathbf{u}}_{\E_0}})$ and the energy of the optimal states $\K({\widetilde{\mathbf{u}}_{\E_0}})$ are shown as functions of $\E_0$ for small $\E_0$ in figures \[fig:RvsE0\_FixE\_small\](a) and \[fig:RvsE0\_FixE\_small\](b), respectively. As indicated by the asymptotic form of $\R$ in and the Poincaré limit [$\K_0=\E_0/(2\pi)^2$]{}, both of which are marked in these figures, the behaviour of $\R_{\E_0}$ and $\K({\widetilde{\mathbf{u}}_{\E_0}})$ as $\E_0 \rightarrow 0$ is correctly captured by the numerically computed optimal states. In particular, we note that $\R_{\E_0}$ is negative for $0 \le
  \E_0 \lessapprox 7$ and exhibits the same trend as predicted in for $\E_0 \rightarrow 0$. For larger values of $\E_0$ the rate of growth of enstrophy becomes positive. Likewise, the asymptotic behaviour of the energy of the optimal fields does not come as a surprise since, as discussed in §\[sec:3D\_InstOpt\_E0to0\], in the limit $\E_0 \to 0$ the maximizers of $\R$ are eigenfunctions of the Laplacian, which also happen to saturate Poincaré’s inequality.

The structure of the optimal vortex states ${\widetilde{\mathbf{u}}_{\E_0}}$ is analyzed next. They are visualized using in which the vortex cores are identified as regions $\Sigma := \{\Gamma_{s}(Q): s
  \ge 0 \}$ for $Q$ defined as [@davidson:turbulence] $$\label{eq:Q_3D}
Q({\mathbf{x}}) := \, \frac{1}{2}\left[ {\operatorname{tr}}(\bm{\Omega}\bm{\Omega}^T) - {\operatorname{tr}}(\mathbf{S}\mathbf{S}^T) \right],$$ where $\mathbf{S}$ and $\bm{\Omega}$ are the symmetric and anti-symmetric parts of the velocity gradient tensor $\nabla{\mathbf{u}}$, that is, $[\mathbf{S}]_{ij} = \frac{1}{2}(\partial_j u_i + \partial_i
u_j)$ and $[\bm{\Omega}]_{ij} = \frac{1}{2}(\partial_j u_i -
\partial_i u_j )$, $i,j=1,2,3$. The quantity $Q$ can be interpreted as the local balance between the strain rate and the vorticity magnitude. The isosurfaces $\Gamma_{0}(Q - 0.5||Q||_{L_\infty})$ representing the optimal states ${\widetilde{\mathbf{u}}_{\E_0}}$ with selected values of $\E_0$ are shown in figures \[fig:RvsE0\_FixE\_small\](c)-(e). For the smallest values of $\E_0$, the optimal fields exhibit a cellular structure already observed in figure \[fig:maxdEdt\_vortexCells\](a). For increasing values of $\E_0$ this cellular structure transforms into a vortex ring, as seen in figure \[fig:RvsE0\_FixE\_small\](e). The component of vorticity normal to the plane $P_x = \{{\mathbf{x}}\in\Omega : \mathbf{n}\cdot({\mathbf{x}}-{\mathbf{x}}_0) = 0
\}$ for $\mathbf{n} = [1,0,0]$ and ${\mathbf{x}}_0 = [1/2,1/2,1/2]$ is shown in figures \[fig:RvsE0\_FixE\_small\](f)-(h), where the transition from cellular structures to a localized vortex structure as enstrophy increases is evident.

The results corresponding to large values of $\E_0$ are shown in figure \[fig:RvsE0\_FixE\_large\] with the maximum rate of growth of enstrophy $\R_{\E_0}$ plotted as a function of $\E_0$ in figure \[fig:RvsE0\_FixE\_large\](a). We observe that, as $\E_0$ increases, this relation approaches a power law of the form $\R_{\E_0} = C_1' \,\E_0^{\alpha_1}$. In order to determine the prefactor $C_1'$ and the exponent $\alpha_1$ we perform a local least-squares fit of the power law to the actual relation $\R_{\E_0}$ versus $\E_0$ for increasing values of $\E_0$ starting with $\E_{0} = 20$ (this particular choice the starting value is justified below). Then, the exponent $\alpha_1$ is computed as the average of the exponents obtained from the local fits with their standard deviation providing the error bars, so that we obtain $$\label{eq:RvsE0_powerLaw}
\R_{\E_0} = C'_1\E_0^{\,\alpha_1}, \qquad C'_1 = 3.72 \times 10^{-3} , \ \alpha_1 = 2.97 \pm 0.02$$ (the same approach is also used to determine the exponents in other power-law relations detected in this study). We note that the exponent $\alpha_1$ obtained in is in fact very close to 3 which is the exponent in estimate . For the value of the viscosity coefficient used in the computations ($\nu=0.01$), the constant factor $C_1 = 27/(8\pi^4\nu^3)$ in estimate has the value $C_1 \approx 3.465 \times 10^4$ which is approximately seven orders of magnitude larger than $C'_1$ given in . To shed more light at the source of this discrepancy, the objective functional $\R$ from equation can be separated into a negative-definite viscous part $\R_{\nu}$ and a cubic part $\R_{\textrm{cub}}$ defined as

$$\begin{aligned}
\R_{\nu}({\mathbf{u}}) & :=  -\nu\int_\Omega |{\Delta}{\mathbf{u}}|^2\,d{\mathbf{x}}, \\
\R_{\textrm{cub}}({\mathbf{u}}) & :=  \int_{\Omega} {\mathbf{u}}\cdot\nabla{\mathbf{u}}\cdot{\Delta}{\mathbf{u}}\, d{\mathbf{x}}, \label{eq:Rcub}\end{aligned}$$

so that $\R({\mathbf{u}}) = \R_{\nu}({\mathbf{u}}) + \R_{\textrm{cub}}({\mathbf{u}})$. The values of $\R_{\textrm{cub}}({\widetilde{\mathbf{u}}_{\E_0}})$ are also plotted in figure \[fig:RvsE0\_FixE\_large\](a) and it is observed that this quantity exhibits a power-law behaviour of the form $$\label{eq:RcubvsE0_powerLaw}
\R_{\textrm{cub}}({\widetilde{\mathbf{u}}_{\E_0}}) = C''_1\E_0^{\,\alpha_2}, \qquad C''_1 =  1.38\times 10^{-2},\ \alpha_2 = 2.99 \pm 0.05.$$ While the value of $C''_1$ is slightly larger than the value of $C'_1$ in , it is still some six orders of magnitude smaller than the constant factor $C_1 = 27/(8\pi^4\nu^3)$ from estimate . These differences notwithstanding, we may conclude that estimate is sharp in the sense of definition \[def:NotionSharpness\]. The power laws from equations and are consistent with the results first presented by @l06 [@ld08], where the authors reported a power-law with exponent $\alpha_{LD} = 2.99$ and a constant of proportionality $C_{LD} = 8.97\times 10^{-4}$. The energy of the optimal fields $\K({\widetilde{\mathbf{u}}_{\E_0}})$ for large values of $\E_0$ is shown in figure \[fig:RvsE0\_FixE\_large\](b) in which we observe that the energy stops to increase at about $\E_0 \approx 20$. This transition justifies using $\E_0 = 20$ as the lower bound on the range of $\E_0$ where the power laws are determined via least-square fits.

Figures \[fig:RvsE0\_FixE\_large\](c)-(e) show the isosurfaces $\Gamma_{0}(Q - 0.5||Q||_{L_\infty})$ representing the optimal fields ${\widetilde{\mathbf{u}}_{\E_0}}$ for selected large values of $\E_0$. The formation of these localized vortex structures featuring two rings as $\E_0$ increases is evident in these figures. The formation process of localized vortex structures is also visible in figures \[fig:RvsE0\_FixE\_large\](f)–(h), where the component of vorticity normal to the plane $P_{xz} = \{{\mathbf{x}}\in\Omega : \mathbf{n}\cdot({\mathbf{x}}- {\mathbf{x}}_0) = 0 \}$ for $\mathbf{n} = [1,0,-1]$ and ${\mathbf{x}}_0 =
[1/2,1/2,1/2]$ is shown (we note that the planes used in figures \[fig:RvsE0\_FixE\_small\](c)–(e) and \[fig:RvsE0\_FixE\_large\](c)–(e) have different orientations).

\
\

\
\

Next we examine the variation of different diagnostics applied to the extreme states ${\widetilde{\mathbf{u}}_{\E_0}}$ as enstrophy $\E_0$ increases. The maximum velocity $||{\widetilde{\mathbf{u}}_{\E_0}}||_{L_\infty}$ and maximum vorticity $||{\bnabla\times}{\widetilde{\mathbf{u}}_{\E_0}}||_{L_\infty}$ of the optimal fields are shown, respectively, in figures \[fig:ScalingLaws\_fixE\](a) and \[fig:ScalingLaws\_fixE\](b) as functions of $\E_0$. For each quantity, two distinct power laws are observed in the forms

$$\begin{aligned}
{4}
||{\widetilde{\mathbf{u}}_{\E_0}}||_{L_\infty} & \sim C_1\E^{\alpha_1}_0,& \quad 
C_1 &= 0.263,&  \alpha_1 & = 0.5 \pm 0.023,&  \quad & 
\mbox{as }\, \E_0\to 0, \\ 
||{\widetilde{\mathbf{u}}_{\E_0}}||_{L_\infty} & \sim C_2\E^{\alpha_2}_0,& \quad 
C_2 &= 6.3\times 10^{-2},& \quad \alpha_2 &= 1.04 \pm 0.13,&  \qquad &
\mbox{as }\, \E_0\to \infty, \label{eq:uLinf}\end{aligned}$$

and

$$\begin{aligned}
{4}
||{\bnabla\times}{\widetilde{\mathbf{u}}_{\E_0}}||_{L_\infty} & \sim C_1\E_0^{\alpha_1},& \quad 
C_1 &= 2.09,& \alpha_1 & = 0.54 \pm 0.03,&  \quad &
\mbox{as }\, \E_0\to 0, \\ 
||{\bnabla\times}{\widetilde{\mathbf{u}}_{\E_0}}||_{L_\infty} & \sim C_2\E_0^{\alpha_2},& \quad 
C_2 &= 6.03\times 10^{-2},& \quad \alpha_2 &= 1.99 \pm 0.17,&  \quad & 
\mbox{as }\, \E_0\to \infty.
\label{eq:omegaLinf}\end{aligned}$$

In order to quantify the variation of the relative size of the vortex structures, we will introduce two characteristic length scales. The first one is based on the energy and enstrophy, and was defined by [@dg95] as $$\label{eq:Lambda_def}
\Lambda := \frac{1}{2\pi}\left[ \frac{\K({\widetilde{\mathbf{u}}_{\E_0}})}{\E({\widetilde{\mathbf{u}}_{\E_0}})} \right]^{1/2}.$$ It is therefore equivalent to the Taylor microscale $\lambda^2
= 15\int_\Omega |{\mathbf{u}}|^2 d{\mathbf{x}}/ \int_\Omega |{\bm{\omega}}|^2 d{\mathbf{x}}$ used in turbulence research [@davidson:turbulence]. Another length scale, better suited to the ring-like vortex structures shown in figures \[fig:RvsE0\_FixE\_large\](c)-(e), is the average radius $R_{\Pi}$ of one of the vortex rings calculated as $$\label{eq:VortexRadius_def}
R_{\Pi} :=  \frac{ \int_\Omega r({\mathbf{x}})\chi_{\Pi}({\mathbf{x}}) \,d{\mathbf{x}}}{ \int_\Omega \chi_{\Pi}({\mathbf{x}})d{\mathbf{x}}},
\ \text{where} \ \
r({\mathbf{x}}) = |{\mathbf{x}}- \overline{{\mathbf{x}}}|,\quad \overline{{\mathbf{x}}} = \frac{\int_{\Omega} {\mathbf{x}}\chi_{\Pi}({\mathbf{x}}) d{\mathbf{x}}}{\int_{\Omega} \chi_{\Pi}({\mathbf{x}})d{\mathbf{x}}},$$ and $\chi_{\Pi}$ is the characteristic function of the set $$\begin{aligned}
\Pi & = & \{ \Gamma_s( Q ) : s > 0.9|| Q ||_{L_\infty}\} \cap \\
 &  & \{ {\mathbf{x}}\in\Omega : {\mathbf{n}}\cdot({\mathbf{x}}-{\mathbf{x}}_0) > 0, \ {\mathbf{n}}= [1,1,1], \ {\mathbf{x}}_0 = [1/2,1/2,1/2] \}.\end{aligned}$$ In the above definition of the set $\Pi$, the intersection of the two regions is necessary to restrict the set $R_{\Pi}$ to only one of the two ring structures visible in figures \[fig:RvsE0\_FixE\_large\](c)–(e). The quantity $\overline{{\mathbf{x}}}$ can be therefore interpreted as the geometric centre of one of the vortex rings. The dependence of $\Lambda$ and $R_\Pi$ on $\E_0$ is shown in figures \[fig:ScalingLaws\_fixE\](c,d) in which the following power laws can be observed

$$\begin{aligned}
{4}
 \Lambda &\sim \O(1)\quad\mbox{and}& \quad  R_\Pi &\sim \O(1)&&& \quad &
 \mbox{as } \,\E_0\to 0, \\
 \Lambda &\sim C_1\E_0^{\alpha_1},& \quad 
 C_1 &= 10.96,&   \alpha_1 &= -0.886 \pm 0.105,& \qquad &
 \mbox{as } \, \E_0\to \infty, \label{eq:Lambda_powerLaw_largeE0} \\ 
 R_\Pi &\sim C_2\E_0^{\alpha_2},& \quad 
 C_2 &= 2.692,& \quad \alpha_2 &= -1.01 \pm 0.16,& \quad &
 \mbox{as } \, \E_0\to \infty. 
\label{eq:Radius_powerLaw_largeE0}\end{aligned}$$

By comparing the error bars in the key power laws and with the error bars in power-law relations , , and , we observe that there is less uncertainty in the first case, indicating that the quantities $||{\widetilde{\mathbf{u}}_{\E_0}}||_{L_\infty}$, $||{\bnabla\times}{\widetilde{\mathbf{u}}_{\E_0}}||_{L_\infty}$, $\Lambda$ and $R_\Pi$ tend to be more sensitive to approximation errors than $\R_{\E_0}({\widetilde{\mathbf{u}}_{\E_0}})$. Non-negligible error bars may also indicate that, due to modest enstrophy values attained in our computations, the ultimate asymptotic regime corresponding to $\E_0 \rightarrow
  \infty$ has not been reached in some power laws.

A useful aspect of employing the average ring radius $R_{\Pi}$ as the characteristic length scale is that its observed scaling with respect to $\E_0$ can be used as an approximate indicator of the resolution $1/N$ required to numerically solve problem \[pb:maxdEdt\_E\] for large values of enstrophy. From the scaling in relation , it is evident that a two-fold increase in the value of $\E_0$ will be accompanied by a similar reduction in $R_\Pi$, thus requiring an eight-fold increase in the resolution (a two-fold increase in each dimension). This is one of the reasons why computation of extreme vortex states ${\widetilde{\mathbf{u}}_{\E_0}}$ for large enstrophy values is a very challenging computational task. In particular, this relation puts a limit on the largest value of $\E_0$ for which problem \[pb:maxdEdt\_E\] can be in principle solved computationally at the present moment: a value of $\E_0 =
  2000$, a mere order of magnitude above the largest value of $\E_0$ reported here, would require a resolution of $8192^3$ used by some of the largest Navier-Stokes simulations to date.

To summarize, as the enstrophy increases from $\E_0 \approx 0$ to $\E_0 = \O(10^2)$, the optimal vortex states change their structure from cellular to ring-like. While with the exception of $\R({\widetilde{\mathbf{u}}_{\E_0}})$ and $\K({\widetilde{\mathbf{u}}_{\E_0}})$, all of the diagnostic quantities behave in a monotonous manner, the corresponding power laws change at about $\E_0
\approx 20$, which approximately marks the transition from the cellular to the ring-like structure (cf. figure \[fig:RvsE0\_FixE\_small\](e) vs. \[fig:RvsE0\_FixE\_large\](c)). This is also the value of the enstrophy beyond which the energy $\K({\widetilde{\mathbf{u}}_{\E_0}})$ starts to decrease (figure \[fig:RvsE0\_FixE\_large\](b)). This transition also coincides with a change of the symmetry properties of the extreme vortex states ${\widetilde{\mathbf{u}}_{\E_0}}$ — while in the limit $\E_0 \rightarrow 0$ these fields feature reflection and discrete rotation symmetries (cf. § \[sec:3D\_InstOpt\_E0to0\]), for $20 \lessapprox \E_0 \rightarrow
  \infty$ the optimal states are characterized by axial symmetry. The asymptotic (as $\E_0 \rightarrow \infty$) extreme vortex states on locally maximizing branches corresponding to the aligned ABC flow and the Taylor-Green vortex (cf. Table \[tab:E0\]) are similar to the fields shown in figures \[fig:RvsE0\_FixE\_large\](c)–(h), except for a different orientation of their symmetry axes with respect to the periodic domain $\Omega$ (these results are not shown here for brevity). The different power laws found here are compared to the corresponding results obtained in 2D in §\[sec:discuss\]. It is also worth mentioning that, as shown by [@ad16], all power laws discussed in this section, cf. , , , and , can be deduced rigorously using arguments based on dimensional analysis under the assumption of axisymmetry for the optimal fields ${\widetilde{\mathbf{u}}_{\E_0}}$.

\

Finally, the findings of this section allow us to shed some light on the “small data” result which provides the conditions on the size of the initial data ${\mathbf{u}}_0$, given in terms of its energy $\K(0)$ and enstrophy $\E(0)$, in the Navier-Stokes system guaranteeing that smooth solutions exist globally in time. The power-law fits and allow us to sharpen condition be replacing the constant on the RHS with either $2\nu / C'_1$ or $2\nu / C''_1$, so that we obtain $$\label{eq:K0E02}
\K(0)\E(0) < \left\{ \frac{2\nu}{C'_1} \ \ \text{or} \ \ \frac{2\nu}{C''_1}  \right\}.$$ The region of the “phase space” $\{\K,\E\}$ described by condition is shown in white in figure \[fig:K0E0\]. The gray region represents the values of $\K(0)$ and $\E(0)$ for which long-time existence of smooth solutions cannot be a priori guaranteed (the two shades of gray correspond to the two constants which can be used in ). Solid circles represent the different extreme states found in this section, whereas the thin curves mark the time-dependent trajectories which will be analyzed in §\[sec:timeEvolution\]. We conclude from figure \[fig:K0E0\] that the change of the properties of the optimal states ${\widetilde{\mathbf{u}}_{\E_0}}$ discussed above occurs in fact for the values of enstrophy ($\E(0) \approx 20$) for which the states ${\widetilde{\mathbf{u}}_{\E_0}}$ are on the boundary of the region of guaranteed long-time regularity.

Time Evolution of Extreme Vortex States {#sec:timeEvolution}
=======================================

The goal of this section is to analyze the time evolution, governed by the Navier-Stokes system , with extreme vortex states identified in §\[sec:3D\_InstOpt\_E\] used as the initial data ${\mathbf{u}}_0$. In particular, we are interested in the finite-time growth of enstrophy $\E(t)$ and how it relates to estimates , and . We will compare these results with the growth of enstrophy obtained using other types of initial data which have also been studied in the context of the blow-up problem for both the Euler and Navier-Stokes systems, namely, the Taylor-Green vortex [@tg37; @bmonmu83; @b91; @bb12], the Kida-Pelz flow [@bp94; @p01; @ghdg08], colliding Lamb-Chaplygin dipoles [@opc12] and perturbed antiparallel vortex tubes [@k93; @k13]. Precise characterization of these different initial conditions is provided in Table \[tab:InitialConditions\] and, for the sake of completeness, the last three states are also visualized in figure \[fig:SummaryIC\]. We comment that, with the exception of the Taylor-Green vortex which was shown in §\[sec:3D\_InstOpt\_E0to0\] to be a local maximizer of problem \[pb:maxdEdt\_E\] in the limit $\E_0 \rightarrow 0$, all these initial conditions were postulated based on rather ad-hoc physical arguments. We also add that, in order to ensure a fair comparison, the different initial conditions listed in Table \[tab:InitialConditions\] are rescaled to have the same enstrophy $\E_0$, which is different from the enstrophy values used in the original studies where these initial conditions were investigated [@opc12; @k13; @dggkpv13; @opmc14]. As regards our choices of the initial enstrophy $\E_0$, to illustrate different possible behaviours, we will consider initial data located in the two distinct regions of the phase space $\{\K,\E\}$ shown in figure \[fig:K0E0\], corresponding to values of $\K_0$ and $\E_0$ for which global regularity may or may not be a priori guaranteed according to estimates –.

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                                                                                            $${\mathbf{u}}_0({\mathbf{x}}) = [u, v, w]$$                                                                                                          
  ------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ---------------------------
                    Instantaneous                                                                                                                                                                                                                                                                                 

     optimizer ${\widetilde{\mathbf{u}}_{\E_0}}$                                                                             $${\mathbf{u}}_0 = \mathop{{\operatorname{argmax}}}_{{\mathbf{u}}\in\mathcal{S}_{\E_0}} \R({\mathbf{u}})$$                                                                                   See problem

                    Taylor-Green                                                                                                                                                                                                                                                                                  

                       vortex                                                                                                                                             $$\begin{aligned}                                                                                                                        $\bm{\gamma}=(1,-1,0)$ in
                                                                                                                                                       u(x,y,z) & = & A\sin(2\pi x)\cos(2\pi y)\cos(2\pi z) \\                                                                                                    
                                                                                                                                                                     v(x,y,z) & = & -A\cos(2\pi x)\sin(2\pi y)\cos(2\pi z) \\                                                                                     
                                                                                                                                                                                         w(x,y,z) & = & 0                                                                                                         
                                                                                                                                                                                 \end{aligned}$$                                                                                                                  

                     equation ,                                                                                                                                                                                                                                                                                   

   $A$ chosen so that $\E({\mathbf{u}}_0) = \E_0$.                                                                                                                                                                                                                                                                

                      Kida-Pelz                                                                                                                                                                                                                                                                                   

                        flow                                                                                                                                              $$\begin{aligned}                                                                                                                               Taken from
                                                                                                                                                     u(x,y,z) & = & A\sin(2\pi x)[ \cos(6\pi y)\cos(2\pi z) - \\                                                                                                  
                                                                                                                                                                                        &   & \cos(2\pi y)\cos(6\pi z) ] \\                                                                                       
                                                                                                                                                                    v(x,y,z) & = & A\sin(2\pi y)[ \cos(6\pi z)\cos(2\pi x) - \\                                                                                   
                                                                                                                                                                                      &   & \cos(2\pi x_3)\cos(6\pi x_1) ] \\                                                                                     
                                                                                                                                                                    w(x,y,z) & = & A\sin(2\pi z)[ \cos(6\pi x)\cos(2\pi y) - \\                                                                                   
                                                                                                                                                                                        &   & \cos(2\pi x)\cos(6\pi y) ] \\                                                                                       
                                                                                                                                                                                 \end{aligned}$$                                                                                                                  

                       @bp94,                                                                                                                                                                                                                                                                                     

   $A$ chosen so that $\E({\mathbf{u}}_0) = \E_0$.                                                                                                                                                                                                                                                                

                   Lamb-Chaplygin                                                                                                                                                                                                                                                                                 

                       dipoles                                                              $$-{\Delta}{\mathbf{u}}_0  =  {\bnabla\times}{\bm{\omega}}_0, \quad {\bm{\omega}}_0 = [\, 0, \, \omega(x,z), \,\omega(x,y)\,]$$ $$\omega(x(r,\theta),y(r,\theta)) = \left\{                                                   Taken from
                                                                                                                                                                            \begin{array}{c@{\,\,}r}                                                                                                              
                                                                                                                                                   -2U\kappa \frac{J_1(\kappa r)}{J_0(\kappa a)}\sin(\theta) & (r \leq a) \\                                                                                      
                                                                                                                                                                                           0 & (r > a)                                                                                                            
                                                                                                                                                                                  \end{array}                                                                                                                     
                                                                                                                                                                                    \right.$$                                                                                                                     

                       @opc12.                                                                                                                                                                                                                                                                                    

            $a = 0.15$, $\kappa a = z_1$,                                                                                                                                                                                                                                                                         

               the first zero of $J_1$                                                                                                                                                                                                                                                                            

        $U = \sqrt{\frac{\E_0}{2\pi z^2_1}}$                                                                                                                                                                                                                                                                      

                      Perturbed                                                                                                                                                                                                                                                                                   

                    anti-parallel                                                                                                                                                                                                                                                                                 

                    vortex tubes                     $$-{\Delta}{\mathbf{u}}_0  =  {\bnabla\times}{\bm{\omega}}_0, \quad {\bm{\omega}}_0 = \omega(x,y)\frac{\bm{\sigma'}}{|\bm{\sigma}'|}(s)$$ $$\omega(x(r,\theta),y(r,\theta)) =  \frac{A}{(r/a)^{4} + 1}$$ $$\bm{\sigma}(s) =  [2a, 2b/\cosh(s^2/c^2)-b, s]$$          Taken from

                        @k13.                                                                                                                                                                                                                                                                                     

                 $a= 0.05$, $b=a/2$,                                                                                                                                                                                                                                                                              

       $c=a$, $s$ is the arc-length parameter.                                                                                                                                                                                                                                                                    

   $A$ chosen so that $\E({\mathbf{u}}_0) = \E_0$.                                                                                                                                                                                                                                                                
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

  : Characterization of the different initial data used in time evolution studies in §\[sec:timeEvolution\].[]{data-label="tab:InitialConditions"}

System is solved numerically with an approach combining a pseudo-spectral approximation of spatial derivatives with a third-order semi-implicit Runge-Kutta method [@NumRenaissance] used to discretize the problem in time. In the evaluation of the nonlinear term dealiasing was used based on the $2/3$ rule together with the Gaussian filtering proposed by [@hl07]. Massively parallel implementation based on MPI and using the [ fftw]{} routines [@fftw] to perform Fourier transforms allowed us to use resolutions varying from $256^3$ to $1024^3$ in the low-enstrophy and high-enstrophy cases, respectively. A number of different diagnostics were checked to ensure that all flows discussed below are well resolved. We refer the reader to the dissertation by [@a14] for additional details and a validation of this approach.

The time-dependent results will be shown with respect to a normalized time defined as $\tau := U_c t/\ell_c$ with $U_c := \|{\widetilde{\mathbf{u}}_{\E_0}}\|_{L_2}$ and $\ell_c = \Lambda$ (cf. equation ) playing the roles of the characteristic velocity and length scale. We begin by showing the time evolution of the enstrophy $\E(\tau)$ corresponding to the five different initial conditions listed in Table \[tab:InitialConditions\] with $\E_0 = 10$ and $\E_0 = 100$ in figures \[fig:Fvs\_t\_fixE\](a) and \[fig:Fvs\_t\_fixE\](b), respectively (because of the faster time-scale, the time axis in the latter figure is scaled logarithmically). We see that the maximizers ${\widetilde{\mathbf{u}}_{\E_0}}$ of problem \[pb:maxdEdt\_E\] are the only initial data which triggers growth of enstrophy for these values of the initial enstrophy and, as expected, this growth is larger for $\E_0 = 100$ than for $\E_0 = 10$. The other initial condition which exhibits some tendency for growth when $\E_0 = 100$ is the Taylor-Green vortex. In all cases the enstrophy eventually decays to zero for large times.

Next we examine whether the flow evolutions starting from the instantaneous maximizers ${\widetilde{\mathbf{u}}_{\E_0}}$ as the initial data saturate the finite-time estimate . We do this by defining functions

\[eq:fg\] $$\begin{aligned}
f(\tau) & := \frac{1}{\E(0)} - \frac{1}{\E(\tau)}\qquad\mbox{and} \label{eq:fga} \\
g(\tau) & := \frac{C}{2\nu}\left[ \K(0) - \K(\tau) \right] \label{eq:fgb}\end{aligned}$$

representing, respectively, the left- and right-hand side of the estimate and then plotting them with respect to the normalized time $\tau$, which is done in figures \[fig:fg\](a) and \[fig:fg\](b) for $\E_0 = 10$ and $\E_0 = 100$, respectively. The constant $C>0$ in the definition of $g(\tau)$ is numerically computed from the power-law fit in . It follows from estimate that $f(\tau) \leq g(\tau)$ pointwise in time. The hypothetical extreme event of a finite-time blow-up can be represented graphically by an intersection of the graph of $f(\tau)$ with the horizontal line $y = 1/\E_0$, which is also shown in figures \[fig:fg\](a)–(b). The behaviour of $g(\tau)$, representing the upper bound in estimate , is quite distinct in figures \[fig:fg\](a) and \[fig:fg\](b) reflecting the fact that the initial data ${\widetilde{\mathbf{u}}_{\E_0}}$ in the two cases comes from different regions of the phase diagram in figure \[fig:K0E0\]. In figure \[fig:fg\](a), corresponding to $\E_0 = 10$, the upper bound $g(\tau)$ never reaches $1/\E_0$, in agreement with the fact that the finite-time blow-up is a priori ruled out in this case. On the other hand, in figure \[fig:fg\](b), corresponding to $\E_0 = 100$, the upper bound $g(\tau)$ does intersect $1/\E_0$ implying that, in principle, finite-time blow-up might be possible in this case. The sharpness of estimate can be assessed by analyzing how closely the behaviour of $f(\tau)$ matches that of $g(\tau)$. In both figures \[fig:fg\](a) and \[fig:fg\](b) we observe that for a short period of time $f(\tau)$ exhibits a very similar growth to the upper bound $g(\tau)$, but then this growth slows down and $f(\tau)$ eventually starts to decrease short of ever approaching the limit $1/\E_0$.

We further characterize the time evolution by showing the maximum enstrophy increase $\delta\E := \mathop{\max}_{t \geq 0} \, \{ \E(t) -
\E(0) \}$ and the time when the maximum is achieved $T_{\max} :=
\mathop{\arg\max}_{t \geq 0}\, \E(t)$ as functions of $\E_0$ in figures \[fig:Emax\_vsE0\_fixE\](a) and \[fig:Emax\_vsE0\_fixE\](b), respectively. In both cases approximate power laws in the form $$\delta\E  \sim  \E^{\alpha_1}_0, \quad \alpha_1 = 0.95 \pm 0.06 \qquad\mbox{and}\qquad 
T_{\max}  \sim  \E^{\alpha_2}_0, \quad \alpha_2  = -2.03 \pm 0.02$$ are detected in the limit $\E_0 \rightarrow \infty$ (as regards the second result, we remark that $T_{\max}$ is not equivalent to the time until which the enstrophy grows at the sustained rate proportional to $\E_0^3$, cf. figure \[fig:fg\]). To complete presentation of the results, the dependence of the quantities $$\mathop{\max}_{t \geq 0} \, \left\{\frac{1}{\E_0} - \frac{1}{\E(t)}\right\} \qquad\mbox{and}\qquad [\K(0) - \K(T_{\max})]$$ on the initial enstrophy $\E_0$ is shown in figures \[fig:Emax\_vsE0\_fixE\](c) and \[fig:Emax\_vsE0\_fixE\](d), respectively. It is observed that both quantities approximately exhibit a power-law behaviour of the form $\E^{-1}_0$. Discussion of these results in the context of the estimates recalled in §\[sec:intro\] is presented in the next section.

\

Discussion {#sec:discuss}
==========

In this section we provide some comments about the results reported in §§\[sec:3D\_InstOpt\_E0to0\], \[sec:3D\_InstOpt\_E\] and \[sec:timeEvolution\]. First, we need to mention that our gradient-based approach to the solution of optimization problem \[pb:maxdEdt\_E\] can only yield local maximizers and, due to nonconvexity of the problem, it is not possible to guarantee a priori that the maximizers found are global. To test for the possible presence of branches other than the ones found using the continuation approach described in §\[sec:3D\_InstOpt\], cf. Algorithm \[alg:optimAlg\], we tried to find new maximizers by initializing the gradient iterations with different initial guesses ${\mathbf{u}}^0$. They were constructed as solenoidal vector fields with prescribed regularity and random structure, which was achieved by defining the Fourier coefficients of ${\mathbf{u}}^0$ as $\widehat{{\mathbf{u}}}^0({\mathbf{k}}) = F(|{\mathbf{k}}|)e^{i\phi({\mathbf{k}})}$ with the amplitude $F(|{\mathbf{k}}|) \sim 1/|{\mathbf{k}}|^2$ and the phases $\phi({\mathbf{k}})$ chosen as random numbers uniformly distributed in $[0,2\pi]$. However, in all such tests conducted for $\E_0 = O(1)$ the gradient optimization algorithm would always converge to maximizers ${\widetilde{\mathbf{u}}_{\E_0}}$ belonging to one of the branches discussed in §\[sec:3D\_InstOpt\_E\] (modulo possible translations in the physical domain). While far from settling this issue definitely, these observations lend some credence to the conjecture that the branch identified in §\[sec:3D\_InstOpt\_E\] corresponds in fact to the global maximizers. These states appear identical to the maximizers found by [@ld08] and our search has also yielded two additional branches of locally maximizing fields, although we did not capture the lower branch reported by [@ld08]. However, since that branch does not appear connected to any state in the limit $\E_0 \rightarrow 0$, we speculate that it might be an artifact of the “discretize-then-optimize” formulation used by [@ld08], in contrast to the “optimize-then-discretize” approach employed in our study which provides a more direct control over the analytic properties of the maximizers. We add that the structure of the maximizing branches found here is in fact quite similar to what was discovered by [@ap13a] in an analogous problem in 2D. Since the 2D problem is more tractable from the computational point of view, in that case we were able to undertake a much more thorough search for other maximizers which did not however yield any solutions not associated with the main branches.

The results reported in §\[sec:3D\_InstOpt\_E\] and §\[sec:timeEvolution\] clearly exhibit two distinct behaviours, depending on whether or not global-in-time regularity can be guaranteed a priori based on estimates –. These differences are manifested, for example, in the power laws evident in figures \[fig:ScalingLaws\_fixE\] and \[fig:Emax\_vsE0\_fixE\], as well as in the different behaviours of the RHS of estimate with respect to time in figures \[fig:fg\](a) and \[fig:fg\](b). However, for the initial data for which global-in-time regularity cannot be ensured a priori there is no evidence of sustained growth of enstrophy strong enough to signal formation of singularity in finite time. Indeed, in figure \[fig:Emax\_vsE0\_fixE\](c) one sees that the quantity $\mathop{\max}_{t \geq 0} \, \left\{1/\E_0 - 1/\E(t)\right\}$ behaves as $C_1 / \E_0$, where $C_1 < 1$, when $\E_0$ increases, revealing no tendency to approach $1/\E_0$ which is a necessary precursor of a singular behaviour (cf. discussion in §\[sec:timeEvolution\]). To further illustrate how the rate of growth of enstrophy achieved initially by the maximizers ${\widetilde{\mathbf{u}}_{\E_0}}$ is depleted in time, in figure \[fig:dEdtE\] we show the flow evolution corresponding to ${\widetilde{\mathbf{u}}_{\E_0}}$ with $\E_0 = 100$ as a trajectory in the coordinates $\{\E, d\E/dt\}$. From the discussion in Introduction we know that in order for the singularity to occur in finite time, the enstrophy must grow at least at a sustained rate $d\E / dt \sim
\E^\alpha$ for some $\alpha > 2$. In other words, a “blow-up trajectory” will be realized only if the trajectory of the flow, expressed in $\{\E,d\E/dt \}$ coordinates, is contained in the region $\mathcal{M} = \{ (\E,d\E/dt) \; : \; C_1\E^2 < d\E/dt \leq C_2\E^3
\}$, for some positive constants $C_1$ and $C_2$. For the flow corresponding to the instantaneous optimizers ${\widetilde{\mathbf{u}}_{\E_0}}$, the initial direction of a trajectory in $\{\E,d\E/dt \}$ coordinates is determined by the vector ${\mathbf{v}}= \left[1,
  \left.\tfrac{d\R}{d\E}\right|_{\E_0} \right]$ and, for initial conditions ${\mathbf{u}}_0$ satisfying $\R({\mathbf{u}}_0) = C\E^3({\mathbf{u}}_0)$, it follows that $$\left.\frac{d\R}{d\E}\right|_{\E_0}  = 3C\E_0^2.$$ Since the optimal rate of growth is sustained only over a short interval of time, the trajectory of the flow in the $\{\E,d\E/dt\}$ coordinates approaches the region $\mathcal{M}$ only tangentially following the direction of the lower bound $C_1\E^2$, and remains outside $\mathcal{M}$ for all subsequent times. This behaviour is clearly seen in the inset of figure \[fig:dEdtE\].

An interpretation of this behaviour can be proposed based on equation from which it is clear that the evolution of the flow energy is closely related to the growth of enstrophy. In particular, if the initial energy $\K(0)$ is not sufficiently large, then its depletion due to the initial growth of enstrophy may render the flow incapable of sustaining this growth over a longer period of time. This is in fact what seems to be happening in the present problem as evidenced by the data shown in figure \[fig:K0E0\]. We remark that, for a prescribed enstrophy $\E(0)$, the flow energy cannot be increased arbitrarily as it is upper-bounded by Poincaré’s inequality [$\K(0) \le (2\pi)^2 \E(0)$]{}. This behaviour can also be understood in terms of the geometry of the extreme vortex states ${\widetilde{\mathbf{u}}_{\E_0}}$. Figure \[fig:ring\] shows a magnification of the pair of vortex rings corresponding to the optimal field ${\widetilde{\mathbf{u}}_{\E_0}}$ with $\E_0=100$. It is observed that the vorticity field ${\bnabla\times}{\widetilde{\mathbf{u}}_{\E_0}}$ inside the vortex core has an azimuthal component only which exhibits no variation in the azimuthal direction. Thus, in the limit $\E_0
\rightarrow \infty$ the vortex ring shrinks with respect to the domain $\Omega$ (cf. figure \[fig:ScalingLaws\_fixE\](d)) and the field ${\widetilde{\mathbf{u}}_{\E_0}}$ ultimately becomes axisymmetric (i.e., in this limit boundary effects vanish). At the same time, it is known that the 3D Navier-Stokes problem on an unbounded domain and with axisymmetric initial data is globally well posed [@k03], a results which is a consequence of the celebrated theorem due to @ckn83.

![Trajectory of the flow corresponding to the initial condition ${\widetilde{\mathbf{u}}_{\E_0}}$ with $\E_0 = 100$ in the coordinates $\{\E,d\E/dt\}$. For comparison, in the inset the thin line represents the borderline growth at the rate $d\E/dt \sim \E^2$](Figs2/dEdt_vsE_NSE_Instant_E100_New){width="90.00000%"}

. \[fig:dEdtE\]

![Vortex lines inside the region with the strongest vorticity in the extreme vortex state ${\widetilde{\mathbf{u}}_{\E_0}}$ with $\E_0 = 100$. The colour coding of the vortex lines is for identification purposes only.[]{data-label="fig:ring"}](Figs2/dnsNS3D_K00_E47_Instant_plane_vorLns){width="60.00000%"}

We close this section by comparing the different power laws characterizing the maximizers ${\widetilde{\mathbf{u}}_{\E_0}}$ and the corresponding flow evolutions with the results obtained in analogous studies of extreme behaviour in 1D and 2D (see also Table \[tab:estimates\]). First, we note that the finite-time growth of enstrophy $\delta\E$ in 3D, cf. figure \[fig:Emax\_vsE0\_fixE\](a), exhibits the same dependence on the enstrophy $\E_0$ of the instantaneously optimal initial data as in 1D, i.e., is directly proportional to $\E_0$ in both cases [@ap11a]. This is also analogous to the maximum growth of palinstrophy $\P$ in 2D which was found by [@ap13a] to scale with the palinstrophy $\P_0$ of the initial data, when the instantaneously optimal initial condition was computed subject to [ *one*]{} constraint only (on $\P_0$). When the instantaneously optimal initial data was determined subject to [*two*]{} constraints, on $\K_0$ and $\P_0$, then the maximum finite-time growth of palinstrophy was found to scale with $\P_0^{3/2}$ [@ap13b]. On the other hand, the time $T_{\max}$ when the maximum enstrophy is attained, cf. figure \[fig:Emax\_vsE0\_fixE\](b), scales as $\E_0^{-2}$, which should be contrasted with the scalings $\E_0^{-1/2}$ and $\P_0^{-1/2}$ found in the 1D and 2D cases, respectively. This implies that the time interval during which the extremal growth of enstrophy is sustained in 3D is shorter than the corresponding intervals in 1D and 2D.

Conclusions and Outlook {#sec:final}
=======================

By constructing the initial data to exhibit the most extreme behaviour allowed for by the mathematically rigorous estimates, this study offers a fundamentally different perspective on the problem of searching for potentially singular solutions from most earlier investigations. Indeed, while the corresponding flow evolutions did not reveal any evidence for finite-time singularity formation, the initial data obtained by maximizing $d\E/dt$ produced a significantly larger growth of enstrophy in finite time than any other candidate initial conditions (cf. Table \[tab:InitialConditions\] and figure \[fig:Fvs\_t\_fixE\]). Admittedly, this observation is limited to the initial data with $\E_0 \le 100$, which corresponds to Reynolds numbers $Re = \sqrt{\E_0\, \Lambda} / \nu \lessapprox 450$ lower than the Reynolds numbers achieved in other studies concerned with the extreme behaviour in the Navier-Stokes flows [@opc12; @k13; @dggkpv13; @opmc14]. Given that the definitions of the Reynolds numbers applicable to the various flow configurations considered in these studies were not equivalent, it is rather difficult to make a precise comparison in terms of specific numerical values, but it is clear that the largest Reynolds numbers attained in these investigations were at least an order of magnitude higher than used in the present study; for Euler flows such a comparison is obviously not possible at all. However, from the mathematical point of view, based on estimates –, there is no clear indication that a very large initial enstrophy $\E_0$ (or, equivalently, a high Reynolds number) should be a necessary condition for singularity formation in finite time. In fact, blow-up cannot be a priori ruled out as soon as condition is violated, which happens for all initial data lying on the gray region of the phase space in figure \[fig:K0E0\]. We remark that additional results were obtained (not reported in this paper) by studying the time evolution corresponding to the optimal initial data ${\widetilde{\mathbf{u}}_{\E_0}}$, but using smaller values of the viscosity coefficient $\nu=10^{-3},
10^{-4}$, thereby artificially increasing the Reynolds number at the price of making the initial data suboptimal. Although these attempts did increase the amplification of enstrophy as compared to what was observed in figures \[fig:fg\] and \[fig:Emax\_vsE0\_fixE\], no signature of finite-time singularity formation could be detected either.

Our study confirmed the earlier findings of [@ld08] about the sharpness of the instantaneous estimate . We also demonstrated that the finite-time estimate is saturated by the flow evolution corresponding to the optimal initial data ${\widetilde{\mathbf{u}}_{\E_0}}$, but only for short times, cf. figure \[fig:fg\], which are not long enough to trigger a singular behaviour.

In §\[sec:discuss\] we speculated that a relatively small initial energy $\K(0)$, cf. figure \[fig:RvsE0\_FixE\_large\](b), might be the property of the extreme vortex states ${\widetilde{\mathbf{u}}_{\E_0}}$ preventing the resulting flow evolutions from sustaining a significant growth of enstrophy over long times. On the other hand, in Introduction we showed that estimate need not be saturated for blow-up to occur in finite time and, in fact, sustained growth at the rate $d\E/dt = C \, \E^\alpha$ with any $\alpha > 2$ will also produce singularity in finite time. Thus, another strategy to construct initial data which could lead to a more sustained growth of enstrophy in finite time might be to increase its kinetic energy by allowing for a smaller instantaneous rate of growth (i.e., with an exponent $2 < \alpha \le 3$ instead of $\alpha = 3$). This can be achieved by prescribing an additional constraint in the formulation of the variational optimization problem, resulting in

\[pb:maxdEdt\_KE\] $$\begin{aligned}
{\widetilde{\mathbf{u}}_{\K_0,\E_0}}& = & \mathop{\arg\max}_{{\mathbf{u}}\in{\mathcal{S}_{\K_0,\E_0}}} \, \R({\mathbf{u}}) \\
{\mathcal{S}_{\K_0,\E_0}} & = & \left\{{\mathbf{u}}\in H_0^2(\Omega)\,\colon\,\nabla\cdot{\mathbf{u}}= 0, \; \K({\mathbf{u}}) = \K_0, \; \E({\mathbf{u}}) = \E_0 \right\}.\end{aligned}$$

It differs from problem \[pb:maxdEdt\_E\] in that the maximizers are sought at the intersection of the original constraint manifold ${\mathcal{S}_{\E_0}}$ and the manifold defined by the condition $\K({\mathbf{u}})
  = \K_0$, where $\K_0 \le (2\pi)^2 \E_0$ is the prescribed energy. While computation of such maximizers is more complicated, robust techniques for the solution of optimization problems of this type have been developed and were successfully used in the 2D setting by [@ap13a]. Preliminary results obtained in the present setting by solving problem \[pb:maxdEdt\_KE\] for $\K_0 =
1$ are indicated in figure \[fig:K0E0\], where we see that the flow evolutions starting from ${\widetilde{\mathbf{u}}_{\K_0,\E_0}}$ do not in fact produce a significant growth of enstrophy either. An alternative, and arguably more flexible approach, is to formulate this problem in terms of multiobjective optimization [@k99] in which the objective function $\R({\mathbf{u}})$ in problem \[pb:maxdEdt\_E\] would be replaced with $$\R_{\eta}({\mathbf{u}}) := \eta\, \R({\mathbf{u}}) + (1-\eta)\, \K({\mathbf{u}}),
\label{eq:Rm}$$ where $\eta \in [0,1]$. Solution of such a multiobjective optimization problem has the form of a “Pareto front” parameterized by $\eta$. Clearly, the limits $\eta \rightarrow 1$ and $\eta \rightarrow 0$ correspond, respectively, to the extreme vortex states already found in §\[sec:3D\_InstOpt\_E0to0\] and §\[sec:3D\_InstOpt\_E\], and to the Poincaré limit. Another interesting possibility is to replace the energy $\K({\mathbf{u}})$ with the helicity [$\H({\mathbf{u}}) :=
  \int_{\Omega} {\mathbf{u}}\cdot(\bnabla\times{\mathbf{u}})\,d\Omega$]{} in the multiobjective formulation , as this might allow one to obtain extreme vortex states with a more complicated topology (i.e., a certain degree of “knottedness”). We note that all the extreme vortex states found in the present study were “unknotted”, i.e., were characterized by $\H({\widetilde{\mathbf{u}}_{\E_0}}) = 0$, as the vortex rings were in all cases disjoint (cf. figure \[fig:ring\]).

Finally, another promising possibility to find initial data producing a larger growth of enstrophy is to solve a [*finite-time*]{} optimization problem of the type already studied by [@ap11a] in the context of the 1D Burgers equation, namely

\[pb:maxdE\] $$\tilde{{\mathbf{u}}}_{0;\E_0,T}  = \mathop{\arg\max}_{{\mathbf{u}}_0\in{\mathcal{S}_{\E_0}}} \, \E(T),$$

where $T>0$ is the length of the time interval of interest and ${\mathbf{u}}_0$ the initial data for the Navier-Stokes system . In contrast to problems \[pb:maxdEdt\_E\] and \[pb:maxdEdt\_KE\], solution of problem \[pb:maxdE\] is more complicated as it involves flow evolution. It represents therefore a formidable computational task for the 3D Navier-Stokes system. However, it does appear within reach given the currently available computational resources and will be studied in the near future.

Acknowledgements {#acknowledgements .unnumbered}
================

The authors are indebted to Charles Doering for many enlightening discussions concerning the research problems studied in this work. The authors are grateful to Nicholas Kevlahan for making his parallel Navier-Stokes solver available, which was used to obtain the results reported in §\[sec:timeEvolution\]. Anonymous referees provided many insightful comments which helped us improve this work. This research was funded through an Early Researcher Award (ERA) and an NSERC Discovery Grant, whereas the computational time was made available by SHARCNET under its Dedicated Resource Allocation Program. Diego Ayala was funded in part by NSF Award DMS-1515161 at the University of Michigan and by the Institute for Pure and Applied Mathematics at UCLA.

[54]{} natexlab\#1[\#1]{}

2005 [*[Sobolev]{} Spaces*]{}. Elsevier.

2014 [Extreme Vortex States and Singularity Formation in Incompressible Flows]{}. PhD thesis, McMaster University, available at .

2016 Extreme enstrophy growth in axisymmetric flows. In preparation.

2016 An estimate for the growth of palinstrophy in [2D]{} incompressible flows. In preparation.

2011 On maximum enstrophy growth in a hydrodynamic system. [*Physica D*]{} [**240**]{}, 1553–1563.

2014 Maximum palinstrophy growth in [2D]{} incompressible flows. [*Journal of Fluid Mechanics*]{} [ **742**]{}, 340–367.

2014 Vortices, maximum growth and the problem of finite-time singularity formation. [*Fluid Dynamics Research*]{} [**46**]{} (3), 031404.

2009 [*Numerical Renaissance*]{}. Renaissance Press.

1994 Direct numerical simulation of transition to turbulence from a high-symmetry initial condition. [*Physics of Fluids*]{} [**6**]{}, 2757–2784.

1991 Direct simulation of three-dimensional turbulence in the [Taylor-Green]{} vortex. [*Fluid Dynamics Research*]{} [**8**]{}, 1–8.

1983 Small-scale structure of the [Taylor-Green]{} vortex. [ *Journal of Fluid Mechanics*]{} [**130**]{}, 411–452.

2012 Interplay between the [Beale-Kato-Majda]{} theorem and the analyticity-strip method to investigate numerically the incompressible [Euler]{} singularity problem. [*Phys. Rev. E*]{} [**86**]{}, 066302.

2008 3d euler about a 2d symmetry plane. [*Physica D*]{} [**237**]{}, 1912–1920.

1983 [Partial regularity of suitable weak solutions of the Navier-Stokes equations]{}. [*Communications on Pure and Applied Mathematics*]{} [**35**]{}, 771–831.

2005 Evolution of complex singularities in [Kida-Pelz]{} and [Taylor-Green]{} inviscid flows. [*Fluid Dynamics Research*]{} [**36**]{}, 239–248.

2004 [*Turbulence. An introduction for scientists and engineers*]{}. Oxford University Press.

2009 The [3D Navier-Stokes]{} problem. [*Annual Review of Fluid Mechanics*]{} pp. 109–128.

1995 [*Applied Analysis of the [Navier-Stokes]{} Equations*]{}. Cambridge University Press.

1986 Chaotic streamlines in the [ABC]{} flows. [*Journal of Fluid Mechanics*]{} [**167**]{}, 353–391.

2013 Vorticity moments in four numerical simulations of the [3D Navier-Stokes]{} equations. [*Journal of Fluid Mechanics*]{} [**732**]{}, 316–331.

2000 Existence and smoothness of the [Navier-Stokes]{} equation. available at , [Clay Millennium Prize Problem Description]{}.

1989 Gevrey class regularity for the solutions of the [Navier–Stokes]{} equations. [*Journal of Functional Analysis*]{} [ **87**]{}, 359–369.

2003 [*[FFTW]{} User’s Manual*]{}. Massachusetts Institute of Technology.

2014 [Regimes of nonlinear depletion and regularity in the 3D Navier-Stokes equations]{}. [*Nonlinearity*]{} [**27**]{} (1–19).

2008 The three–dimensional [Euler]{} equations: singular or non–singular? [*Nonlinearity*]{} [**21**]{}, 123–129.

2008 Numerical simulations of possible finite-time singularities in the incompressible [Euler]{} equations: comparison of numerical methods. [*Physica D*]{} [ **237**]{}, 1932–1936.

2003 [*Perspectives in Flow Control and Optimization*]{}. SIAM.

2009 Blow-up or no blow-up? a unified computational and analytic approach to [3D]{} incompressible [Euler]{} and [Navier–Stokes]{} equations. [*Acta Numerica*]{} pp. 277–346.

2007 Computing nearly singular solutions using pseudo-spectral methods. [*Journal of Computational Physics*]{} [**226**]{}, 379–397.

1993 Evidence for a singularity of the three-dimensional, incompressible [Euler]{} equations. [*Phys. Fluids A*]{} [**5**]{}, 1725–1746.

2013 Bounds for [Euler]{} from vorticity moments and line divergence. [*Journal of Fluid Mechanics*]{} [**729**]{}, R2.

2013 Swirling, turbulent vortex rings formed from a chain reaction of reconnection events. [*Physics of Fluids*]{} [**25**]{}, 065101.

2003 [Remarks for the axisymmetric Navier-Stokes equations]{}. [ *Journal of Differential Equations*]{} [**187**]{}, 226–239.

2010 Regularity and blow up for active scalars. [*Math. Model. Nat. Phenom.*]{} [**5**]{}, 225–255.

2004 [*Initial-Boundary Value Problems and the [Navier-Stokes]{} Equations*]{}, [*Classics in Applied Mathematics*]{}, vol. 47. SIAM.

1969 [*The mathematical Theory of Viscous Incompressible Flow*]{}. Gordon and Breach.

2006 Bounds on the enstrophy growth rate for solutions of the [3D Navier-Stokes]{} equations. PhD thesis, University of Michigan.

2008 Limits on enstrophy growth for solutions of the three-dimensional [Navier–Stokes]{} equations. [*Indiana University Mathematics Journal*]{} [**57**]{}, 2693–2727.

1969 [*Optimization by Vector Space Methods*]{}. John Wiley and Sons.

2014 [Potentially Singular Solutions of the 3D Axisymmetric Euler Equations]{}. [*Proceedings of the National Academy of Sciences*]{} (in press).

2014 [Toward the Finite-Time Blowup of the 3D Incompressible Euler Equations: a Numerical Investigation]{}. [*SIAM: Multiscale Modeling and Simulation*]{} (in press).

2002 [*Vorticity and Incompressible Flow*]{}. Cambridge University Press.

2008 Complex-space singularities of [2D Euler]{} flow in lagrangian coordinates. [*Physica D*]{} [**237**]{}, 1951–1955.

1999 [*Nonlinear Multiobjective Optimization*]{}. Springer.

2008 A miscellany of basic issues on incompressible fluid equations. [*Nonlinearity*]{} [**21**]{}, 255–271.

2008 [Numerical study of the Eulerian–Lagrangian analysis of the [Navier-Stokes]{} turbulence]{}. [*Phys. Fluids*]{} [**20**]{}, 1–11.

2014 A minimal flow unit for the study of turbulence with passive scalars. [ *Journal of Turbulence*]{} [**15**]{}, 731–751.

2012 Vortex events in [Euler]{} and [Navier-Stokes]{} simulations with smooth initial conditions. [ *Journal of Fluid Mechanics*]{} [**690**]{}, 288–320.

2001 Symmetry and the hydrodynamic blow-up problem. [ *Journal of Fluid Mechanics*]{} [**444**]{}, 299–320.

2004 A comprehensive framework for the regularization of adjoint analysis in multiscale [PDE]{} systems. [ *Journal of Computational Physics*]{} [**195**]{}, 49–89.

1990 Collapsing solutions to the [3D Euler]{} equations. [*Phys. Fluids A*]{} [**2**]{}, 220–241.

2006 [*Nonlinear Optimization*]{}. Princeton University Press.

2009 Calculation of complex singular solutions to the [3D]{} incompressible [Euler]{} equations. [*Physica D*]{} [**238**]{}, 2368–2379.

1937 Mechanism of the production of small eddies from large ones. [*Proceedings of the Royal Society of London A*]{} [**158**]{}, 499–521.

[^1]: Email address for correspondence: bprotas@mcmaster.ca
