#  Systematic search for singularities in 3D Euler flows

Xinyu Zhao1 and Bartosz Protas1, Corresponding Author, Email:
<EMAIL_ADDRESS>
( 1Department of Mathematics and Statistics, McMaster University
Hamilton, ON, Canada
)

###### Abstract

We consider the question whether starting from a smooth initial condition 3D
inviscid Euler flows on a periodic domain $\mathbb{T}^{3}$ may develop
singularities in a finite time. Our point of departure is the well-known
result by Kato (1972), which asserts the local existence of classical
solutions to the Euler system in the Sobolev space $H^{m}(\mathbb{T}^{3})$ for
$m>5/2$. Thus, potential formation of a singularity must be accompanied by an
unbounded growth of the $H^{m}$ norm of the velocity field as the singularity
time is approached. We perform a systematic search for “extreme” Euler flows
that may realize such a scenario by formulating and solving a PDE-constrained
optimization problem where the $H^{3}$ norm of the solution at a certain fixed
time $T>0$ is maximized with respect to the initial data subject to suitable
normalization constraints. This problem is solved using a state-of-the-art
Riemannian conjugate gradient method where the gradient is obtained from
solutions of an adjoint system. Computations performed with increasing
numerical resolutions demonstrate that, as asserted by the theorem of Kato
(1972), when the optimization time window $[0,T]$ is sufficiently short, the
$H^{3}$ norm remains bounded in the extreme flows found by solving the
optimization problem, which indicates that the Euler system is well-posed on
this “short” time interval. On the other hand, when the window $[0,T]$ is
long, possibly longer than the time of the local existence asserted by Kato’s
theorem, then the $H^{3}$ norm of the extreme flows diverges upon resolution
refinement, which indicates a possible singularity formulation on this “long”
time interval. The extreme flow obtained on the long time window has the form
of two colliding vortex rings and is characterized by certain symmetries. In
particular, the region of the flow in which a singularity might occur is
nearly axisymmetric.

##  1 Introduction

The motion of an ideal (inviscid) incompressible fluid on a domain
$\Omega\subseteq\mathbb{R}^{3}$ is described by the Euler equations

$\displaystyle\frac{\partial\bm{u}}{\partial
t}+(\bm{u}\cdot\bm{\nabla})\bm{u}$ $\displaystyle=-\nabla p,$
$\displaystyle(\bm{x},t)\in\Omega\times(0,T],$ (1.1)
$\displaystyle\bm{\nabla}\cdot\bm{u}$ $\displaystyle=0,$
$\displaystyle(\bm{x},t)\in\Omega\times(0,T],$ $\displaystyle\bm{u}|_{t=0}$
$\displaystyle=\bm{\eta},$ $\displaystyle\bm{x}\in\Omega,$

where $\bm{\nabla}=[\partial_{1},\partial_{2},\partial_{3}]^{T}$ is the
gradient operator, $\bm{u}=\bm{u}(\bm{x},t)=[u_{1},u_{2},u_{3}]^{T}$ is the
velocity field, $p=p(\bm{x},t)$ is the scalar pressure and $T>0$ is the length
of the time window considered. We assume the flow domain to be a periodic unit
cube $\Omega=\mathbb{T}^{3}:=\mathbb{R}^{3}/\mathbb{Z}^{3}$, where “:=” means
“equal to by definition”. In system (1.1), the symbol $\bm{\eta}$ denotes the
initial condition for the velocity which is assumed to satisfy
$\bm{\nabla}\cdot\bm{\eta}=0$. We will use the notation $\bm{u}(t;\bm{\eta})$
to indicate the dependence of the solution of system (1.1) at time $t$ on
$\bm{\eta}$. We are interested in the possibility of a spontaneous formation
of finite-time singularities in solutions of the Euler system (1.1) equipped
with smooth initial data $\bm{\eta}$, the so-called “blow-up problem”. It
remains one of the central open questions in mathematical fluid mechanics
(Gibbon et al., 2008) and is closely related to the corresponding regularity
problem for the viscous Navier-Stokes system (Doering, 2009; Robinson, 2020),
which has been recognized by the Clay Mathematics Institute as one of its
“millennium problems” (Fefferman, 2000). Our goal in this paper is to
undertake a systematic search for possible singularities in the Euler system
(1.1) using methods of numerical PDE optimization.

The first results asserting local existence of classical solutions to the
Euler system were obtained by Lichtenstein (1925) for Hölder-regular initial
data $\bm{\eta}\in C^{1,\alpha}(\mathbb{R}^{3})$ with $0<\alpha<1$. Local
existence of classical solutions in Sobolev spaces was then established by
Kato (1972) and is summarized in the following theorem

###### Theorem 1.1

If $\bm{\eta}\in H^{m}(\mathbb{T}^{3})$ for some $m>5/2$ and satisfies
$\bm{\nabla}\cdot\bm{\eta}=0$, then there exists a time
$T=T\left(\|\bm{\eta}\|_{H^{m}}\right)>0$ such that (1.1) has a unique
solution $\bm{u}(\cdot;\bm{\eta})\in C([0,T];H^{m})\bigcap
C^{1}([0,T];H^{m-1})$.

The study of the local well-posedness of the Euler system (1.1) with analytic
initial data began with the work of Bardos & Benachour (1977).

Another well-known conditional regularity result is the Beale-Kato-Majda (BKM)
criterion (Beale et al., 1984) which states that a smooth solution $\bm{u}$ of
the Euler system develops a singularity at $t=T^{*}$ if and only if

$\int_{0}^{T^{*}}||\bm{\omega}(t)||_{L^{\infty}}\;dt=\infty,$ (1.2)

where $\bm{\omega}:=\nabla\times\bm{u}$ is the vorticity field of the flow.
The sufficiency of this condition can be deduced from Theorem 1.1 using a
Sobolev inequality (Adams & Fournier, 2005)

$\int_{0}^{T^{*}}||\bm{\omega}(t)||_{L^{\infty}}\;dt\leq T^{*}\sup_{0\leq
t\leq T^{*}}||\nabla\bm{u}(t)||_{L^{\infty}}\leq CT^{*}\sup_{0\leq t\leq
T^{*}}||\bm{u}(t)||_{H^{m}},\quad m>5/2,$ (1.3)

whereas its necessity is a result of the inequality (Beale et al., 1984; Majda
& Bertozzi, 2002)

$||\bm{u}(t)||_{H^{m}}\leq||\bm{\eta}||_{H^{m}}\exp\left\\{C_{1}\exp\left(C_{2}\int_{0}^{T^{*}}||\bm{\omega}(\tau)||_{L^{\infty}}\;d\tau\right)\right\\}.$
(1.4)

The double exponential in the upper bound in (1.4) suggests that, should blow-
up indeed occur at some $t=T^{*}$, we can expect a much rapid growth of
$\|\bm{u}(t)\|_{H^{m}}$, $m>5/2$, than of $\|\bm{\omega}(t)\|_{L^{\infty}}$,
as $t\rightarrow T^{*}$. There have been various refinements of the BKM
criterion, where the $L^{\infty}$ norm of vorticity in (1.2) is replaced with
a norm in the BMO space (Kozono & Taniuchi, 2000) or a Besov space (Chae,
2001). In addition, there are also geometric criteria for blow-up (Constantin
et al., 1996; Deng et al., 2005), in which the direction of the vorticity
$\bm{\omega}/|\bm{\omega}|$ plays a crucial role. As was shown in Gibbon
(2013), the relative ordering of suitably rescaled vorticity moments provides
information about the degree of depletion of the nonlinearity, and hence also
about the regularity of solutions. On the other hand, regularity of weak
solutions of the Euler system is related to Onsager’s conjecture concerning
energy dissipation in such flows (Constantin et al., 1994) and significant
progress has been made recently as regards the non-uniqueness of weak
solutions with low regularity (De Lellis & Székelyhidi Jr, 2019).

With respect to blow-up scenarios, Elgindi (2021) proved that there exist
swirl-free axisymmetric solutions of the Euler system corresponding to initial
data in $C^{1,\alpha}(\mathbb{R}^{3})$ that can develop finite-time
singularities. Elgindi & Jeong (2019) showed the existence of finite-time
singularities in strong solutions to axisymmetric Euler equations on exterior
domains with “hourglass”-shaped boundaries. Most recently, Chen & Hou (2022)
proved a nearly self-similar blow-up of solutions to the 3D axisymmetric Euler
equations with smooth initial data on cylindrical domains with solid
boundaries. We emphasize that in the last two cases the presence of the solid
boundary is key to the formation of the singularity. A simplified semi-
analytic model of vortex reconnection and blow-up was recently developed and
analyzed based on the Biot-Savart law and asymptotic techniques by Moffatt &
Kimura (2019a, b).

The computational studies exploring the possibility of finite-time blow-up in
the Euler system include Brachet et al. (1983); Pumir & Siggia (1990); Brachet
(1991); Kerr (1993); Pelz (2001); Bustamante & Kerr (2008); Ohkitani &
Constantin (2008); Ohkitani (2008); Grafke et al. (2008); Gibbon et al.
(2008); Hou (2009); Orlandi et al. (2012); Bustamante & Brachet (2012); Kerr
(2013); Orlandi et al. (2014); Campolina & Mailybaev (2018), all of which
considered problems subject to periodic boundary conditions in all three
spatial directions. We also mention the studies by Matsumoto et al. (2008) and
Siegel & Caflisch (2009), along with references found therein, in which
various complexified forms of the Euler equations were investigated. The idea
of this approach is that, since 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. Some of the investigations (Bustamante & Brachet, 2012; Orlandi
et al., 2012) hinted at the possibility of singularity formation a in finite
time. In this connection we also highlight the computational investigations of
Luo & Hou (2014a, b) in which blow-up was documented in axisymmetric Euler
flows on a bounded cylindrical domain. This mechanism of singularity formation
involves an interaction with the solid boundary and was recently validated
with rigorous mathematical analysis by Chen & Hou (2022). We also mention an
investigation by Hou (2022) who provided evidence for blow-up in axisymmetric
Euler flows on bounded domains in which singularity occurs away from the
boundaries. In contrast to most other studies, the works of Luo & Hou (2014a,
b); Hou (2022) relied on adaptive mesh refinement employed to resolve fine
structures in flows at the edge of regularity. In this context, we also
mention the investigation by Yin et al. (2021) where the Euler system (1.1)
was solved in the Lagrangian setting using the characteristic mapping method.
Regarding weak solutions and Onsager’s related conjecture, Fehn et al. (2022)
recently presented numerical evidence for nonvanishing energy dissipation in
weak solutions of the Euler system (1.1) which was obtained using a method
based on a discontinuous Galerkin approximation.

While in the aforementioned computational studies the initial conditions for
the Euler system were chosen in an ad-hoc manner, albeit one usually motivated
by deep physical considerations, here we follow a fundamentally different
approach wherein the initial data leading to the most singular, in a
mathematically precise sense, solutions is sought systematically by solving a
suitably constrained PDE optimization problem. This approach was originally
proposed by Lu & Doering (2008) as a way to determine incompressible velocity
fields maximizing the instantaneous rate of growth of enstrophy in 3D Navier-
Stokes flows. It was later extended by Ayala & Protas (2011, 2014, 2017), Yun
& Protas (2018) to study extreme behavior in solutions of different
hydrodynamic models such as variants of the Burgers equation. Recently, Kang
et al. (2020) and Kang & Protas (2021) applied this method to search for
potential singularities in 3D Navier-Stokes flows based on some classical
conditional regularity results. Highlights of this more than decade-long
research program are summarized in the review paper by Protas (2022). The
present investigation represents a first application of this framework to an
inviscid problem.

Guided by the local well-posedness result in Theorem 1.1, we aim to find an
initial condition $\bm{\eta}\in H^{m}$ subject to certain constraints, such
that the $H^{m}$ norm of the corresponding solution of the Euler system (1.1)
is maximized at a prescribed time $T$. The desired initial conditions are thus
found as local maximizers of a constrained PDE optimization problem with the
square of the $H^{m}$ (semi)norm used as the objective functional where for
concreteness we set $m=3$. We solve this optimization problem for different
time intervals $[0,T]$ using a Riemannian conjugate gradient method (Absil et
al., 2008), where the gradient is conveniently computed from solutions of a
suitably-defined adjoint system. Both the Euler system (1.1) and the adjoint
system as well as the different diagnostic quantities are approximated
numerically using pseudospectral methods.

In analogy to the work of Fehn et al. (2022); Guo et al. (2022), we adopt an
indirect approach to distinguish between regular and singular evolution based
on resolution refinement. When solving the optimization problem on a “short”
time interval $[0,T]$, if the objective functional approximated using
different resolutions converges to a finite value as the resolution is
refined, then we conclude the Euler system (1.1) is well-posed on this short
interval. However, when the interval $[0,T]$ is “long”, presumably longer than
the minimum time of existence guaranteed by Theorem 1.1, the objective
functional evaluated on the optimal solutions may diverge upon resolution
refinements. In practice, based on numerical experiments, we use $T=25$ and
$T=75$ for the short and long optimization time windows. We do in fact observe
these two distinct behaviors and conclude that in the latter scenario the
divergence of the objective functional in Euler flows on the long time
interval $[0,75]$ may signal a potential singularity formation in the flow.
The behavior of the different diagnostic quantities does not contradict the
possibility of a singularity formation, although a firm conclusion cannot be
drawn due to insufficient numerical resolution. The corresponding flow
features two colliding vortex rings and while these flow structures are quite
deformed, the flow does exhibit certain nontrivial symmetries. In particular,
the region in the flow where the potential singularity may occur is nearly
axisymmetric.

The structure of the paper is as follows: in § 2 we formulate an optimization
problem designed to eludicate the most extreme behavior possible in an Euler
flow and in § 3 we describe a Riemannian gradient-based approach we use to
solve this problem together with its numerical discretization; the extreme
flows found in this way on a “short” and “long” time intervals are discussed
in § 4; finally, conclusions and outlook are deferred to § 5.

##  2 Optimization problem

Before introducing the optimization problem, let us first define the function
space in which its solutions will be sought. We begin by defining the norm in
the Sobolev space $H^{m}(\mathbb{T}^{3})$ as (Adams & Fournier, 2005)

$\|\bm{u}\|_{H^{m}}:=\left[\sum_{\bm{j}\in\mathbb{Z}^{3}}\left(1+|2\pi\bm{j}|^{2}\right)^{m}|\hat{\bm{u}}_{\bm{j}}|^{2}\right]^{1/2},\qquad
m\in\mathbb{R},$ (2.1)

where $\hat{\bm{u}}_{\bm{j}}$ is the Fourier coefficient of $\bm{u}$
corresponding to the wavevector $\bm{k}:=2\pi\bm{j}=[2\pi j_{1},2\pi
j_{2},2\pi j_{3}]$. The corresponding seminorm, denoted
$\|\cdot\|_{\dot{H}^{m}}$, is obtained from (2.1) by dropping the constant
($\bm{j}$-independent) terms.

###  2.1 Functional setting

Since system (1.1) is locally well posed in $H^{m}(\Omega)$, $m>5/2$, cf.
Theorem 1.1, it may appear natural to look for optimal initial data
$\bm{\eta}$ in that space. However, we are interested in finite-time
singularities potentially arising in classical solutions, whereas initial
conditions constructed in such spaces will in general not be smooth or real-
analytic. Moreover, solving the Euler system (1.1) with such initial data
would not allow us to benefit from the exponential convergence of the
pseudospectral methods used in this study, cf. § 3.3.1. We will thus consider
an extended Gevrey space $G^{\sigma}$ with $\sigma>0$ of real-analytic
functions defined on $\mathbb{T}^{3}$ and will endow it with the inner product

$\displaystyle\forall\bm{v},\bm{u}\in
G^{\sigma},\quad\left\langle\bm{v},\bm{u}\right\rangle_{G^{\sigma}}:=$
$\displaystyle\sum_{\bm{j}\in\mathbb{Z}^{3}}(1+|2\pi\bm{j}|^{2})^{m}e^{4\pi\sigma|\bm{j}|}\hat{\bm{v}}_{\bm{j}}\cdot\overline{\hat{\bm{u}}}_{\bm{j}}$
(2.2) $\displaystyle=$
$\displaystyle\int_{\mathbb{T}^{3}}(1+|D|^{2})^{m}e^{2\sigma|D|}\bm{v}\cdot\bm{u}\,d\bm{x},$

where overbar denotes complex conjugation, whereas the operators $|D|$ and
$e^{\sigma|D|}$ are defined via

$\left[\widehat{|D|\bm{v}}\right]_{\bm{j}}:=2\pi|\bm{j}|\hat{\bm{v}}_{\bm{j}},\qquad\qquad\left[\widehat{e^{\sigma|D|}\bm{v}}\right]_{\bm{j}}:=e^{2\pi\sigma|\bm{j}|}\hat{\bm{v}}_{\bm{j}}.$
(2.3)

In the setting of our problem, the Gevrey space can be regarded a linear
subspace of the Sobolev space $H^{m}(\Omega)$, i.e.,

$G^{\sigma}:=\left\\{\bm{v}\in H^{m}(\Omega)\ :\
\|\bm{v}\|_{G^{\sigma}}=\left\langle\bm{v},\bm{v}\right\rangle^{1/2}_{G^{\sigma}}<\infty\right\\},\quad
m>\frac{5}{2},\quad\sigma>0.$ (2.4)

###### Remark 2.1

It follows from definitions (2.2) and (2.4) that the Gevrey spaces have the
property

$G^{\sigma_{2}}\subset G^{\sigma_{1}}\subset G^{0}=H^{m},\qquad
0<\sigma_{1}<\sigma_{2}.$ (2.5)

Our strategy for choosing the value of $\sigma$ will be discussed in § 4.1.

Since the initial condition $\bm{\eta}$ in system (1.1) needs to be
divergence-free and the quantity $\int_{\mathbb{T}}\bm{\eta}\,d\bm{x}$ is an
invariant of motion, we introduce the following subspace in which optimal
initial conditions will be sought

$V:=\\{\bm{v}\in
G^{\sigma}:\,\bm{\nabla}\cdot\bm{v}=0,\quad\int_{\mathbb{T}^{3}}\bm{v}\,d\bm{x}=\bm{0}\\}.$
(2.6)

The Gevrey regularity is closely related to the analyticity of complex
extensions of functions in $G^{\sigma}$, i.e.,
$\bm{u}(\bm{z}):=\bm{u}(\bm{x}+i\bm{y})$, where $i:=\sqrt{-1}$. If $\bm{u}\in
G^{\sigma}$ can be extended to a strip
$S_{\delta}:=\\{\bm{z}=[z_{1},z_{2},z_{3}]^{T}\in\mathbb{C}^{3}:|\operatorname{Im}(\bm{z}_{1})|+|\operatorname{Im}(\bm{z}_{2})|+|\operatorname{Im}(\bm{z}_{3})|<\delta\\}$,
so that $\bm{u}(\bm{z})$ does not have any singularities in $S_{\delta}$, then
the largest such value of $\delta$ is called the width of the analyticity
strip and we have $0<\sigma\leq\delta$. Furthermore, if $\bm{u}=\bm{u}(t)$ is
a smooth time-dependent solution of the Euler system (1.1), then the width of
the analyticity strip is a function of time, $\delta=\delta(t)$.

As regards analyticity of solutions of the Euler system (1.1), among related
works we mention the study by Bardos & Benachour (1977) which showed that if
the initial data $\bm{\eta}$ is analytic, the solution $\bm{u}(t;\bm{\eta})$
also remains analytic as long as it is well defined in a certain Hölder space.
The relation between the BKM condition (1.2) and the Gevrey regularity was
investigated by Kukavica & Vicol (2009) who obtained lower bounds on
$\delta(t)$ proportional to
$\exp\left(-\int_{0}^{t}||\nabla\bm{u}(\cdot,s)||_{L^{\infty}}\;ds\right)$.
The width of the analyticity strip $\delta(t)$ can be conveniently estimated
from the Fourier spectrum of the solution $\bm{u}(t)$ (Sulem et al., 1983) and
is often used as a regularity indicator whose vanishing implies a singularity
formation in numerical computations of Euler flows (Bustamante & Brachet,
2012).

###  2.2 Statement of the optimization problem

Motivated by the local existence result in Theorem 1.1, the goal of our
optimization-based formulation is to construct initial data $\bm{\eta}$ with a
prescribed $H^{m}$, $m>5/2$, norm, such that at a given time $T>0$, the
corresponding solution $\bm{u}(T;\bm{\eta})$ of the Euler system (1.1) has the
largest possible $H^{m}$ norm. To fix attention, we will hereafter use $m=3$.
Furthermore, since we consider zero-mean initial conditions, cf. (2.6),
without loss of generality, we can use the $H^{3}$ seminorm, denoted
$\|\cdot\|_{\dot{H}^{3}}$, instead of the full $H^{3}$ norm as this will make
the optimization problem more similar to the problem considered by Kang et al.
(2020).

We define the objective functional $\Phi_{T}:V\to\mathbb{R}^{+}$ as

$\Phi_{T}(\bm{\eta}):=\|\bm{u}(T;\bm{\eta})\|^{2}_{\dot{H}^{3}},$ (2.7)

where $\bm{u}(T;\bm{\eta})$ is the solution of the Euler system (1.1) obtained
with the initial condition $\bm{\eta}\in V$. System (1.1) possesses a scaling
property such that if $\left\\{\bm{u}(\bm{x},t),p(\bm{x},t)\right\\}$ is a
solution, then for any $\lambda>0$,

$\left\\{\bm{u}^{\lambda}(\bm{x},t):=\lambda\bm{u}(\bm{x},\lambda t),\qquad
p^{\lambda}(\bm{x},t):=\lambda^{2}p(\bm{x},\lambda t)\right\\}$ (2.8)

is also a solution. Hence, for any nonzero initial condition $\bm{u}_{0}$, we
have the identity

$\Phi_{T}(\bm{\eta})=||\bm{\eta}||^{2}_{\dot{H}^{3}}\Phi_{||\bm{\eta}||_{\dot{H}^{3}}T}\left(\frac{\bm{\eta}}{||\bm{\eta}||_{\dot{H}^{3}}}\right).$
(2.9)

Therefore, we can restrict our discussion to initial conditions with unit
$\dot{H}^{3}$ seminorm which belong to a closed manifold
$\mathcal{M}_{1}\subset V$ defined as

$\mathcal{M}_{1}=\\{\bm{v}\in V:||\bm{v}||_{\dot{H}^{3}}=1\\}.$ (2.10)

Thus, we arrive at the following optimization problem

###### Problem 2.2

Given $T\in\mathbb{R}_{+}$, find

$\widetilde{\bm{\eta}}_{T}=\operatorname*{argmax}\limits_{\bm{\eta}\in\mathcal{M}_{1}}\Phi_{T}(\bm{\eta}).$
(2.11)

Flows corresponding to optimal initial conditions $\widetilde{\bm{\eta}}_{T}$
are referred to as “extreme”. Since a priori we do not know if a singularity
may appear in Euler flows on a given time interval $[0,T]$, we solve the
optimization problem for increasing values of $T$ and deduce whether a
singularity may occur on the interval $[0,T]$ by performing resolution
refinement as discussed in detail in § 4. If the time window $[0,T]$ falls
within the interval of local existence guaranteed by Theorem 1.1, then for any
$\bm{\eta}\in\mathcal{M}_{1}$, we expect $\Phi_{T}(\bm{\eta})$ to be finite
such that it will remain bounded upon resolution refinement. On the other
hand, if there exists an $\bm{\eta}\in\mathcal{M}_{1}$, which will lead to a
finite-time blow-up inside a sufficiently long time interval $[0,T]$, we
anticipate $\Phi_{T}(\bm{\eta})$ to diverge as the resolution is refined.

##  3 Solution approach

We solve Problem 2.2 using a Riemannian conjugate gradient approach which is a
modification of a gradient ascent method given by the iterative relation

$\bm{\eta}_{T}^{(n+1)}=\bm{\eta}_{T}^{(n)}+\tau_{n}\nabla\Phi_{T}\left(\bm{\eta}_{T}^{(n)}\right),\qquad\bm{\eta}_{T}^{(0)}=\bm{\eta}_{0},\qquad
n=0,1,\dots,$ (3.1)

where $\nabla\Phi_{T}\left(\bm{\eta}^{(n)}\right)$ is the gradient of the
objective functional $\Phi_{T}(\bm{\eta})$ evaluated at the element
$\bm{\eta}_{T}^{(n)}$, $\tau_{n}$ is the step size along the gradient
direction and $\bm{\eta}_{0}$ is the initial guess. A local maximizer
$\widetilde{\bm{\eta}}_{T}$ of Problem 2.2 can then be found as the limit
$\widetilde{\bm{\eta}}_{T}=\lim_{n\to\infty}\bm{\eta}_{T}^{(n)}$, such that

$\Phi_{T}\left(\widetilde{\bm{\eta}}_{T}\right)=\lim_{n\to\infty}\Phi_{T}\left(\bm{\eta}_{T}^{(n)}\right)=\max_{\bm{\eta}\in\mathcal{M}_{1}}\Phi_{T}(\bm{\eta})=:\widetilde{\Phi}_{T;\;\bm{\eta}_{0}}.$
(3.2)

Since our optimization problem is non-convex, the sequence constructed in
(3.1) may converge to different local maximizers depending on the initial
guess $\bm{\eta}_{0}$. Therefore, we use the subscript $\bm{\eta}_{0}$ on the
right-hand side (RHS) in (3.2) to indicate the dependence of the local
maximizer on the initial guess.

In the spirit of the “optimize-then-discrtize” paradigm (Gunzburger, 2003), we
first formulate our approach in the infinite-dimensional (continuous) setting
and then discretize the resulting relations in numerical computations. A key
element of the iterative procedure (3.1) is the evaluation of the gradient
$\nabla\Phi_{T}\left(\bm{\eta}\right)$, which is based on solving a suitably-
defined adjoint system backwards in time as discussed in § 3.1. The Riemannian
conjugate gradient method, which accelerates gradient ascent method (3.1) for
problems defined on smooth manifolds such as $\mathcal{M}_{1}$, is introduced
in § 3.2, whereas numerical techniques used to discretize the problem are
briefly discussed in § 3.3.1. Finally, different initial guesses we use are
listed in § 3.3.2 and in § 3.3.3 we describe the diagnostic quantities that
will be analyzed.

###  3.1 Evaluation of the gradient

To compute the gradient $\nabla\Phi_{T}(\bm{\eta})$ of the objective
functional in (3.1), we first introduce the Gâteaux (directional) differential
$\Phi_{T}^{\prime}(\bm{\eta};\bm{\eta}^{\prime}):V\times V\to\mathbb{R}$

$\Phi_{T}^{\prime}(\bm{\eta};\bm{\eta}^{\prime}):=\lim_{\epsilon\to
0}\frac{1}{\epsilon}\left[\Phi_{T}\left(\bm{\eta}+\epsilon\bm{\eta}^{\prime}\right)-\Phi_{T}(\bm{\eta})\right],$
(3.3)

which represents the variation of the objective functional
$\Phi_{T}(\bm{\eta})$ resulting from applying an infinitesimal perturbation
proportional to $\bm{\eta}^{\prime}\in V$ to the initial condition
$\bm{\eta}$. Fixing the first argument of
$\Phi_{T}^{\prime}(\bm{\eta};\bm{\eta}^{\prime})$, we can view the Gâteaux
differential (3.3) as a bounded linear functional on $V$. Therefore, we can
define $\nabla\Phi_{T}(\bm{\eta})$ using the Riesz representation theorem
(Luenberger, 1969) in terms of the inner product on $G^{\sigma}$, cf. (2.4),
as

$\Phi_{T}^{\prime}(\bm{\eta};\bm{\eta}^{\prime})=\left\langle\nabla\Phi_{T}(\bm{\eta}),\bm{\eta}^{\prime}\right\rangle_{G^{\sigma}},\qquad\bm{\eta},\bm{\eta}^{\prime}\in
V.$ (3.4)

The Gâteaux differential can be evaluated by substituting (2.7) into (3.3),
which yields

$\Phi_{T}^{\prime}(\bm{\eta};\bm{\eta}^{\prime})=2\left\langle\bm{u}(\cdot,T),\bm{u}^{\prime}(\cdot,T)\right\rangle_{\dot{H}^{3}}=2\left\langle|D|^{6}\bm{u}(\cdot,T),\bm{u}^{\prime}(\cdot,T)\right\rangle_{L^{2}}.$
(3.5)

Here $\bm{u}^{\prime}(\bm{x},t)$ is the solution of the linearization of the
Euler system (1.1) around its solution $\bm{u}(\bm{x},t;\bm{\eta})$, which is
defined by

$\displaystyle\mathcal{L}\begin{bmatrix}\bm{u}^{\prime}\\\
p^{\prime}\end{bmatrix}:$
$\displaystyle=\begin{bmatrix}\partial_{t}\bm{u}^{\prime}+\bm{u}^{\prime}\cdot\bm{\nabla}\bm{u}+\bm{u}\cdot\bm{\nabla}\bm{u}^{\prime}+\bm{\nabla}p^{\prime}\\\\[3.0pt]
\bm{\nabla}\cdot\bm{u}^{\prime}\end{bmatrix}=\begin{bmatrix}\bm{0}\\\\[3.0pt]
0\end{bmatrix},$ (3.6) $\displaystyle\bm{u}^{\prime}(0)$
$\displaystyle=\bm{\eta}^{\prime},$

where $\bm{\eta}^{\prime}$ is the perturbation of the initial condition and
$p^{\prime}(\bm{x},t)$ is the corresponding perturbation pressure.

We note, however, that expression (3.5) for the Gâteaux differential is not
yet consistent with the Riesz form (3.4), since the perturbation
$\bm{\eta}^{\prime}$ of the initial condition does not appear in it
explicitly, but is instead “hidden” in the initial condition of the
perturbation problem (3.6). In order to transform expression (3.5) into the
required Riesz form, where the perturbation $\bm{\eta}^{\prime}$ appears
explicitly as a linear factor, we introduce the adjoint states
$\left\\{\bm{u}^{*}:\mathbb{T}^{3}\times[0,T]\to\mathbb{R}^{3},\
p^{*}:\mathbb{T}^{3}\times[0,T]\to\mathbb{R}\right\\}$ and the following
duality-pairing relation

$\displaystyle\left(\mathcal{L}\begin{bmatrix}\bm{u}^{\prime}\\\
p^{\prime}\end{bmatrix},\,\begin{bmatrix}\bm{u}^{*}\\\
p^{*}\end{bmatrix}\right):=$
$\displaystyle\int_{0}^{T}\int_{\mathbb{T}^{3}}\mathcal{L}\begin{bmatrix}\bm{u}^{\prime}\\\
p^{\prime}\end{bmatrix}\cdot\begin{bmatrix}\bm{u}^{*}\\\
p^{*}\end{bmatrix}\,d\bm{x}\,dt$ (3.7) $\displaystyle=$
$\displaystyle\left(\begin{bmatrix}\bm{u}^{\prime}\\\
p^{\prime}\end{bmatrix},\,\mathcal{L}^{*}\begin{bmatrix}\bm{u}^{*}\\\
p^{*}\end{bmatrix}\right)+\int_{\mathbb{T}^{3}}\bm{u}^{\prime}(\bm{x},T)\cdot\bm{u}^{*}(\bm{x},T)\,d\bm{x}$
$\displaystyle-\int_{\mathbb{T}^{3}}\bm{\eta}^{\prime}(\bm{x})\cdot\bm{u}^{*}(\bm{x},0)\,d\bm{x}$
$\displaystyle=$ $\displaystyle 0.$

Here “$\cdot$” denotes the usual Euclidean inner product on $\mathbb{R}^{4}$
and $\mathcal{L}^{*}$ is the adjoint operator defined in terms of the
following system, which is obtained by performing integration by parts with
respect to both space and time in (3.7)

$\displaystyle\mathcal{L}^{*}\begin{bmatrix}\bm{u}^{*}\\\ p^{*}\end{bmatrix}:$
$\displaystyle=\begin{bmatrix}-\partial_{t}\bm{u}^{*}-(\bm{\nabla}\bm{u}^{*}+(\bm{\nabla}\bm{u}^{*})^{T})\bm{u}-\bm{\nabla}p^{*}\\\\[3.0pt]
-\bm{\nabla}\cdot\bm{u}^{*}\end{bmatrix}=\begin{bmatrix}\bm{0}\\\\[3.0pt]
0\end{bmatrix},$ (3.8) $\displaystyle\bm{u}^{*}(T)$
$\displaystyle=2|D|^{6}\bm{u}(T).$

We remark that system (3.8) is a terminal-value problem and as such has to be
integrated backwards in time. We also note that the spatial mean of its
solutions is conserved during the time evolution and therefore we have

$\int_{\mathbb{T}^{3}}\bm{u}^{*}(\bm{x},T)\;d\bm{x}=\int_{\mathbb{T}^{3}}\bm{u}^{*}(\bm{x},0)\;d\bm{x}=0.$
(3.9)

Combining (3.4), (3.5), (3.7) and (3.8), we have

$\Phi^{\prime}_{T}(\bm{\eta};\bm{u}^{\prime}_{0})=\int_{\mathbb{T}^{3}}\bm{u}^{\prime}(\bm{x},T)\cdot\bm{u}^{*}(\bm{x},T)\,d\bm{x}=\int_{\mathbb{T}^{3}}\bm{u}^{\prime}_{0}(\bm{x})\cdot\bm{u}^{*}(\bm{x},0)\,d\bm{x}=\left\langle\nabla\Phi_{T}(\bm{\eta}),\bm{\eta}^{\prime}\right\rangle_{G^{\sigma}}.$
(3.10)

Using the definition of the Gevrey inner product (2.2) in the last equality
above, we obtain

$\displaystyle\left\langle\nabla\Phi_{T}(\bm{\eta}),\bm{\eta}^{\prime}\right\rangle_{G^{\sigma}}$
$\displaystyle=\int_{\mathbb{T}^{3}}\bm{\eta}^{\prime}(\bm{x})\cdot\bm{u}^{*}(\bm{x},0)\,d\bm{x}$
(3.11)
$\displaystyle=\int_{\mathbb{T}^{3}}(1+|D|^{2})^{3}e^{2\sigma|D|}\left\\{\bm{\eta}^{\prime}(\bm{x})\cdot(1+|D|^{2})^{-3}e^{-2\sigma|D|}\bm{u}^{*}(\bm{x},0)\right\\}\;d\bm{x}$
$\displaystyle=\left\langle(1+|D|^{2})^{-3}e^{-2\sigma|D|}\bm{u}^{*}(\bm{x},0),\bm{\eta}^{\prime}(\bm{x})\right\rangle_{G^{\sigma}},$

which allows us to identify the gradient as

$\nabla\Phi_{T}(\bm{\eta})=(1+|D|^{2})^{-3}e^{-2\sigma|D|}\bm{u}^{*}(0).$
(3.12)

In summary, evaluation of this expression requires the solution of the adjoint
system (3.8), which depends on the solution of the Euler system (1.1) with the
initial condition $\bm{\eta}$.

###  3.2 Riemannian conjugate gradient method

For any $\bm{\eta}\in\mathcal{M}_{1}$,
$\bm{\eta}+\tau\nabla\Phi_{T}(\bm{\eta})$ given by (3.1) is divergence-free
but does not necessarily have a unit $\dot{H}^{3}$ seminorm and thus will not
in general belong to $\mathcal{M}_{1}$. In order to enforce this constraint
and accelerate iterations, we use a three-step procedure: first, we project
the gradient $\nabla\Phi_{T}(\bm{\eta})$ onto the subspace
$\mathcal{T}_{\bm{\eta}}\mathcal{M}_{1}$ tangent to $\mathcal{M}_{1}$ at
$\bm{\eta}$, then perform a suitable vector transport operation in order to
construct a Riemannian conjugate ascent direction using the previous search
direction and, finally, retract the resulting state back to the constraint
manifold $\mathcal{M}_{1}$.

We first construct the projection onto a subspace
$\mathcal{T}_{\bm{v}}\mathcal{M}_{1}$ tangent to the constraint manifold at
$\bm{v}\in\mathcal{M}_{1}$ which is defined as

$\mathcal{T}_{\bm{v}}\mathcal{M}_{1}=\Big{\\{}\bm{z}\in
G^{\sigma}(\mathbb{T}^{3}):\,\bm{\nabla}\cdot\bm{z}=0,\quad\int_{\mathbb{T}^{3}}\bm{z}\,d\bm{x}=\bm{0},\quad\langle\bm{v},\bm{z}\rangle_{\dot{H}^{3}}=0\Big{\\}},$
(3.13)

where the last condition represents the Gâteaux differential of the expression
$\|\bm{v}\|_{\dot{H}^{3}(\Omega)}=1$. For any
$z\in\mathcal{T}_{\bm{v}}\mathcal{M}_{1}$, this condition can be expressed as

$\displaystyle\langle\bm{v},\bm{z}\rangle_{\dot{H}^{3}}$
$\displaystyle=\int_{\mathbb{T}^{3}}|D|^{3}\bm{v}\cdot|D|^{3}\bm{z}\,d\bm{x}$
(3.14)
$\displaystyle=\int_{\mathbb{T}^{3}}(1+|D|^{2})^{3}e^{2\sigma|D|}\Big{(}(1+|D|^{2})^{-3}e^{-2\sigma|D|}|D|^{6}\bm{v}\Big{)}\cdot\bm{z}\,d\bm{x}$
$\displaystyle=\left\langle(1+|D|^{2})^{-3}e^{-2\sigma|D|}|D|^{6}\bm{v},\bm{z}\right\rangle_{G^{\sigma}}=0.$

Thus, at each point $\bm{v}\in\mathcal{M}_{1}$, the tangent space
$\mathcal{T}_{\bm{v}}\mathcal{M}_{1}$ can be characterized using the unit
“normal” vector given by

$\bm{n}_{\bm{v}}=\frac{(1+|D|^{2})^{-3}e^{-2\sigma|D|}|D|^{6}\bm{v}}{\big{|}\big{|}(1+|D|^{2})^{-3}e^{-2\sigma|D|}|D|^{6}\bm{v}\big{|}\big{|}_{G^{\sigma}}}$
(3.15)

which allows us to construct the projection operator
$\operatorname{P}_{\mathcal{T}_{\bm{v}}\mathcal{M}_{1}}\;:\;G^{\sigma}\rightarrow{\mathcal{T}_{\bm{v}}\mathcal{M}_{1}}$,
cf. Figure 1,

$\forall\bm{w}\in
G^{\sigma},\qquad\operatorname{P}_{\mathcal{T}_{\bm{v}}\mathcal{M}_{1}}\bm{w}=\bm{w}-\left\langle\bm{n}_{\bm{v}},\bm{w}\right\rangle_{G^{\sigma}}\bm{n}_{\bm{v}}.$
(3.16)

In addition, we also need to introduce the retraction operator
$\operatorname{R}\;:\;{\mathcal{T}_{\bm{v}}\mathcal{M}_{1}}\rightarrow\mathcal{M}_{1}$
which has the form of normalization needed to satisfy the constraint
$||\bm{v}||_{\dot{H}^{3}}=1$ (Absil et al., 2008; Sato, 2021)

$\operatorname{R}(\bm{v})=\frac{\bm{v}}{||\bm{v}||_{\dot{H}^{3}}},\qquad\bm{v}\neq\bm{0}.$
(3.17)

In the standard (Euclidean) conjugate gradient approach, the search direction
at $(n+1)$-th iteration is obtained as a suitable linear combination of the
steepest ascent direction (i.e., the gradient) computed at $(n+1)$-th
iteration and the search direction obtained at $n$-th iteration. In the
Riemannian setting the difficulty with this approach is that these two
directions do not belong to the same linear space and hence cannot be added
directly. In order to circumvent this difficulty, we introduce the vector
transport operation defined as

${{\cal TM}_{1}\oplus{\cal TM}_{1}\rightarrow{\cal
TM}_{1}\;:\;(\bm{\varphi},\bm{\xi})\longmapsto\Gamma_{\bm{\varphi}}(\bm{\xi})\in{\cal
TM}_{1},}$ (3.18)

where ${\cal TM}_{1}=\cup_{\bm{v}\in{\cal M}_{1}}{{\cal T}_{\bm{v}}{\cal
M}_{1}}$ is the tangent bundle on $\mathcal{M}_{1}$, describing how the vector
field $\bm{\xi}$ is transported along the manifold $\mathcal{M}_{1}$ in the
direction of $\bm{\varphi}$ (Absil et al., 2008). The vector transport
provides a map between the tangent spaces ${{\cal
T}_{\bm{\eta}_{T}^{(n)}}{\cal M}_{1}}$ and ${{\cal
T}_{\bm{\eta}_{T}^{(n+1)}}{\cal M}_{1}}$ obtained at two consecutive
iterations, so that algebraic operations can be performed on vectors belonging
to these subspaces. It therefore generalizes the concept of the parallel
translation to the motion on the manifold. In general, vector transport is
defined up to a multiplicative prefactor and can be computed via the
differentiated retraction as (Absil et al., 2008)

$\mathcal{T}_{\bm{\varphi}_{\bm{u}}}(\bm{\xi}_{\bm{u}})=\frac{d}{dt}\operatorname{R}_{\bm{u}}({\bm{\varphi}_{\bm{u}}+t\bm{\xi}_{\bm{u}}})\big{|}_{t=0}=\frac{1}{\|\bm{u}+\bm{\varphi}_{\bm{u}}\|}\left[\bm{\xi}_{\bm{u}}-\frac{\langle\bm{u}+\bm{\varphi}_{\bm{u}},\bm{\xi}_{\bm{u}}\rangle}{\|\bm{u}+\bm{\varphi}_{\bm{u}}\|^{2}}(\bm{u}+\bm{\varphi}_{\bm{u}})\right],\qquad\bm{u}\in\mathcal{M}_{1}$
(3.19)

Figure 1: Schematic illustration of the Riemannian conjugate gradient method
(3.20).

Combining the elements introduced above, the Riemannian conjugate gradient
approach is given by the iterative relation

$\bm{\eta}_{T}^{(n+1)}=\operatorname{R}\Big{[}\bm{\eta}_{T}^{(n)}+\tau_{n}\bm{d}^{(n)}\Big{]},\qquad\bm{\eta}_{T}^{(0)}=\bm{\eta}_{0},\qquad
n=0,1,\dots,$ (3.20)

and is also schematically illustrated in Figure 1. Here the search direction
$\bm{d}^{(n)}$ is computed as

$\displaystyle\bm{d}^{(0)}$
$\displaystyle=\operatorname{P}_{\mathcal{T}_{0}}\bm{G}^{(0)},$ (3.21)
$\displaystyle\bm{d}^{(n)}$
$\displaystyle=\operatorname{P}_{\mathcal{T}_{n}}\bm{G}^{(n)}+\beta_{n}\Gamma_{\tau_{n-1}\bm{d}^{n-1}}\left(\bm{d}^{n-1}\right),\qquad
n\geq 1,$

where $G^{(n)}:=\nabla\Phi_{T}\left(\bm{\eta}^{(n)}\right)$ and
$\mathcal{T}_{n}:=\mathcal{T}_{\bm{\eta}^{(n)}}$. The “momentum” term
$\beta_{n}$ is chosen to enforce the conjugacy of consecutive search
directions and is computed using the Polak-Ribière approach (Nocedal & Wright,
1999)

$\beta_{n}=\frac{\Big{\langle}\operatorname{P}_{\mathcal{T}_{n}}\bm{G}^{(n)},\Big{(}\operatorname{P}_{\mathcal{T}_{n}}\bm{G}^{(n)}-\Gamma_{\tau_{n-1}\bm{d}^{n-1}}\operatorname{P}_{\mathcal{T}_{n-1}}\bm{G}^{(n-1)}\Big{)}\Big{\rangle}_{G^{\sigma}}}{\Big{|}\Big{|}\operatorname{P}_{\mathcal{T}_{n-1}}\bm{G}^{(n-1)}\Big{|}\Big{|}^{2}_{G^{\sigma}}}.$
(3.22)

The step size $\tau_{n}$ in (3.20) is determined by solving the arc-search
problem

$\tau_{n}=\operatorname*{argmax}_{\tau>0}\left\\{\Phi_{T}\left(\operatorname{R}\Big{[}\bm{\eta}^{(n)}+\tau\bm{d}^{(n)}\Big{]}\right)\right\\}$
(3.23)

using a suitable derivative-free approach, such as a variant of Brent’s
algorithm (Press et al., 1986). Due to the presence of the retraction operator
(3.17) in (3.20), the search in (3.23) is performed following a geodesic arc
on the manifold $\mathcal{M}_{1}$, this problem can be regarded as a
generalization of the more common line search approach (Press et al., 1986).
An application of a Riemannian conjugate gradient method to solve an infinite-
dimensional minimization problem with a similar structure arising in quantum
fluids is described by (Danaila & Protas, 2017).

###  3.3 Numerical methods

####  3.3.1 Discretization, validation and numerical parameters

In numerical computations, systems (1.1) and (3.8) are discretized in space
using a standard pseudospectral Fourier method where derivatives are evaluated
in the Fourier space, whereas nonlinear products are computed in the physical
space on an uniform grid with $N$ points in each direction. As regards
dealiasing, we use both the 3/2 rule and the Gaussian filter proposed by Hou &
Li (2007). The resulting system of ordinary differential equations is
discretized in time using an explicit fourth-order Runge-Kutta (RK4) method.
The code is parallelized using MPI and the fast Fourier transforms (FFTs) are
computed using the parallel version of the software library FFTW (Frigo &
Johnson, 2003). We consider four spatial resolutions with
$N^{3}=128^{3},256^{3},512^{3},1024^{3}$ and the time steps $\Delta t$ are
adjusted to satisfy the stability conditions of the RK4 method. If needed, the
resolution used will be denoted with a superscript $N$, i.e., $\bm{u}^{N}(t)$
represents the numerical solution of the Euler system (1.1) computed with the
resolution $N^{3}$, etc. At each iteration in (3.20), the solution $\bm{u}(t)$
of the Euler system (1.1) is needed to construct coefficients in the adjoint
system (3.8). For the three smaller resolutions, we store the solutions of
(1.1) at discrete time steps so that they can be read directly from files when
solving (3.8). However, this is not possible for the resolution of $1024^{3}$
due to storage limitations. We get around this difficulty using the time-
reversibility of the Euler system (1.1) to evolve this system backwards in
time while solving the adjoint system; this allows us to determine the
coefficients in the latter system without having to store the entire solution
to the former.

The solver for the Euler system (1.1) was validated by verifying the
convergence of the numerical solutions to certain exact solutions of (1.1) as
we increase the spatial resolution and/or decrease the time step. As regards
the accuracy of the evaluation of the gradient
$\nabla\Phi_{T}\left(\bm{\eta}\right)$, we consider the quantity (Ayala, 2014)

$\kappa^{N}(\epsilon):=\frac{\Phi_{T}^{N}\left(\bm{\eta}+\epsilon\bm{\eta}^{\prime}\right)-\Phi_{T}^{N}(\bm{\eta})}{\epsilon\left\langle\nabla\Phi_{T}^{N}(\bm{\eta}),\bm{\eta}^{\prime}\right\rangle_{G^{\sigma}}}$
(3.24)

which is the ratio of a finite-difference approximation of the Gâteaux
differential (3.3) and its representation in terms of the Riesz formula (3.4),
both approximated using the resolution $N^{3}$. In theory, we expect that
$\lim_{N\rightarrow\infty}\lim_{\epsilon\rightarrow 0}\kappa^{N}(\epsilon)=1$.
However, in numerical computations performed with finite-precision
arithmetics, $\kappa^{N}(\epsilon)$ will become unbounded as
$\epsilon\rightarrow 0$ due to subtractive cancellation errors in the
numerator in (3.24). For intermediate values of $\epsilon$ the quantity
$|\kappa^{N}(\epsilon)-1|$ is a measure of inaccuracies involved in the
computation of $\Phi_{T}^{N}\left(\bm{\eta}\right)$ and
$\nabla\Phi_{T}^{N}\left(\bm{\eta}\right)$ due to approximation errors in the
numerical solution of systems (1.1) and (3.8). In our validation tests we
observed that, as expected, $|\kappa^{N}(\epsilon)-1|$ is reduced as the
discretization parameters $N$ and $\Delta t$ are refined.

In computations with the nonlinear conjugate gradient method, it is common to
restart the algorithm periodically with a simple gradient step, i.e., by
setting $\beta_{n}=0$, which allows us to improve convergence by getting rid
of ineffectual old information (Nocedal & Wright, 1999). In our computations
we restart the algorithm (3.20) based on the following two criteria:

  * (1)

$n=20k$, $k\in\mathbb{Z}_{+}$,

  * (2)

The search direction $\bm{d}^{(n)}$ fails to be an ascent direction, i.e.,

$\frac{\Big{\langle}\bm{d}^{(n)},\operatorname{P}_{\mathcal{T}_{n}}G^{(n)}\Big{\rangle}_{G^{\sigma}}}{\Big{|}\Big{|}\bm{d}^{(n)}\Big{|}\Big{|}_{G^{\sigma}}\Big{|}\Big{|}\operatorname{P}_{\mathcal{T}_{n}}G^{(n)}\Big{|}\Big{|}_{G^{\sigma}}}<\epsilon_{\text{machine}},\qquad
0<\epsilon_{\text{machine}}\ll 1.$ (3.25)

Iterations (3.20) are declared converged when the relative change of the
objective functional between two consecutive iterations becomes smaller than a
specified tolerance $0<\text{tol}\ll 1$, i.e., when

$\frac{\widetilde{\Phi}_{T}^{N}\left(\bm{\eta}_{T}^{(n+1)}\right)-\widetilde{\Phi}_{T}^{N}\left(\bm{\eta}_{T}^{(n)}\right)}{\widetilde{\Phi}_{T}^{N}\left(\bm{\eta}_{T}^{(n)}\right)}<\text{tol}.$
(3.26)

####  3.3.2 Initial guesses

Since Problem 2.2 is non-convex, iterations (3.20) may produce local
maximizers which will depend on the initial guess $\bm{\eta}_{0}$. In
numerical computations, we thus consider a variety of different initial
guesses listed below, all normalized such that
$\|\bm{\eta}_{0}\|_{\dot{H}^{3}(\Omega)}=1$. For simplicity, we use the same
notation for the original initial conditions and their normalizations. In each
case, we solve Problem 2.2 to obtain an optimal initial condition
$\widetilde{\bm{\eta}}_{T}$ such that for any given $T$ the norm
$\|\bm{u}(T;\widetilde{\bm{\eta}}_{T})\|_{\dot{H}^{3}}$ is larger than
$\|\bm{u}(T;\bm{\eta}_{0})\|_{\dot{H}^{3}}$ where $\bm{\eta}_{0}$ is each of
the following initial guesses.

  * (1)

The 3D Taylor-Green vortex (Taylor & Green, 1937):

$\bm{\eta}_{{\text{TG}}}:=\begin{bmatrix}\sin(2\pi x)\cos(2\pi y)\cos(2\pi
z)\\\ -\cos(2\pi x)\sin(2\pi y)\cos(2\pi z)\\\ 0\end{bmatrix}.$ (3.27)

The Taylor-Green vortex has been widely used as a potential candidate for
blow-up in Euler flows, however, it is still an open question whether this
initial condition can indeed lead to a singularity formation in a finite time
(Cichowlas & Brachet, 2005; Bustamante & Brachet, 2012; Fehn et al., 2022). We
mention that the time windows $[0,T]$ considered here are much shorter than
the times when a potentially singular behavior was observed in these earlier
studies.

  * (2)

Random initial condition:

${\bm{\eta}_{{\text{rand}}}}:=\operatorname{P}_{\text{L}}{\bm{v}_{\text{rand}}}:={\bm{v}_{\text{rand}}}-\bm{\nabla}\Delta^{-1}(\bm{\nabla}\cdot{\bm{v}_{\text{rand}}}),$
(3.28)

where $\operatorname{P}_{\text{L}}$ is the Leray projector (Majda & Bertozzi,
2002) with the property that for any $\bm{w}\in H^{1}(\Omega)$,
$\bm{\nabla}\cdot\left(\operatorname{P}_{\text{L}}\bm{w}\right)=0$ and
$\bm{v}_{\text{rand}}$ is given by

${\bm{v}_{\text{rand}}}=\begin{bmatrix}v_{1}\\\ v_{2}\\\
v_{3}\end{bmatrix}=\begin{bmatrix}\sum_{|j_{1}|+|j_{2}|+|j_{3}|\leq
N_{0}}e^{-|\bm{j}|/(4\pi)}e^{i2\pi\theta_{1}(j_{1})}e^{i2\pi\bm{j}\cdot\bm{x}}+c.c.\\\\[3.0pt]
\sum_{|j_{1}|+|j_{2}|+|j_{3}|\leq
N_{0}}e^{-|\bm{j}|/(4\pi)}e^{i2\pi\theta_{2}(j_{2})}e^{i2\pi\bm{j}\cdot\bm{x}}+c.c.\\\\[3.0pt]
\sum_{|j_{1}|+|j_{2}|+|j_{3}|\leq
N_{0}}e^{-|\bm{j}|/(4\pi)}e^{i2\pi\theta_{3}(j_{3})}e^{i2\pi\bm{j}\cdot\bm{x}}+c.c.\end{bmatrix},$
(3.29)

in which $\theta_{1}$, $\theta_{2}$ and $\theta_{3}$ are $N_{0}$-dimensional
random variables with uniform distributions on $[0,1]^{N_{0}}$. In practice,
we choose $N_{0}=64$ and emphasize that in contrast to the unimodal Taylor-
Green vortex (3.27), the random initial condition (3.28) has energy
distributed over Fourier coefficients with wavevectors within the shell of
radius $2\pi N_{0}$.

  * (3)

Kerr’s initial condition (Kerr, 1993; Hou & Li, 2006):

${\bm{\eta}_{\text{K}}}:=\bm{\nabla}\times(|D|^{-2}\bm{\omega}_{{\text{K}}})$
(3.30)

which represents two perturbed anti-parallel vortex tubes located
symmetrically with respect to the $xy$-plane such that
$\bm{\omega}_{{\text{K}}}(x,y,z)=-\bm{\omega}_{{\text{K}}}(x,y,-z)$. Here the
vorticity field is rescaled to the domain $[-1/2,1/2]^{3}$ and the vorticity
of the vortex tube above the $xy$-plane is given by

$\bm{\omega}_{{\text{K}}}=8G\omega(r)[\omega_{1},\omega_{2},\omega_{3}]^{T},$
(3.31)

where $G$ is a Fourier filter defined by
$\hat{G}_{\bm{j}}=\exp\left(-0.05\left(j_{1}^{4}+j_{2}^{4}+j_{3}^{4}\right)\right)$
and

$\begin{gathered}{\omega(r)=\begin{cases}\exp\left[\frac{-r^{2}}{1-r^{2}}+r^{4}(1+r^{2}+r^{4})\right],&r<1,\\\
0,&r\geq 1,\end{cases}}\quad
r=\frac{1}{R}\sqrt{[x-x(s)]^{2}+[z-z(s)]^{2}},\\\\[3.0pt]
x(s)=\delta_{x}\cos(\pi s/L_{x}),\quad z(s)=z_{0},\\\\[3.0pt]
s(y)=y_{2}+L_{y}\delta_{y_{1}}\sin(\pi y_{2}/L_{y}),\quad y_{2}(y)=4\pi
y+L_{y}\delta_{y_{2}}\sin(\pi(4\pi y)/L_{y}),\\\\[3.0pt]
\omega_{1}=-\frac{\pi\delta_{x}}{L_{x}}\left[1+\pi\delta_{y_{2}}\cos(\pi(4\pi
y)/L_{y})\right]\cdot\left[1+\pi\delta_{y_{1}}\cos(\pi
y_{2}/L_{y})\right]\sin(\pi s/L_{x}),\\\\[3.0pt]
\omega_{2}=1,\qquad\qquad\omega_{3}=0.\end{gathered}$ (3.32)

We choose the same parameters as Hou & Li (2006) and Kerr (1993), i.e.,

$\delta_{y_{1}}=0.5,\;\delta_{y_{2}}=0.4,\;\delta_{x}=-1.6,\;z_{0}=1.57,\;R=0.75,\;L_{x}=L_{y}=4\pi,\;L_{z}=2\pi.$
(3.33)

We would like to point out that in Hou & Li (2006), the equations are posed on
$[-2\pi,2\pi]^{3}$ rather than $[-1/2,1/2]^{3}$. We also note that with the
normalization used here the times up to which Kerr (1993); Hou & Li (2006)
were able to compute solutions of the Euler system (1.1) with the initial
condition $\bm{\eta}_{0}=\bm{\eta}_{\text{K}}$ corresponding to, respectively,
$t_{1}=108231$ and $t_{2}=110059$, which are significantly longer than the
optimization windows $T$ we consider in this study.

  * (4)

Hou’s axisymmetric initial condition (Hou, 2022):

$\bm{\eta}_{\text{H}}:=Gu_{\theta}\bm{e}_{\theta},$ (3.34)

where $\bm{e}_{\theta}$ is the unit vector in the azimuthal direction of the
cylindrical coordinate system and $u_{\theta}$ is the angular velocity
component given by

$u_{\theta}=\begin{cases}r\exp\left(-r^{2}/(1-r^{2})\right)\frac{12000(1-r^{2})^{18}\sin(2\pi
z)}{1+12.5\sin^{2}(\pi z)},&r<1,\\\ 0,&r\geq 1,\end{cases}\quad
r=\frac{\sqrt{x^{2}+y^{2}}}{0.9}.$ (3.35)

This initial condition is a “compactified” version of the initial condition
originally considered by Hou (2022), where it was defined on a cylindrical
domain of unit radius and with a no-flow boundary condition. The results
reported by Hou (2022) indicate that the Euler flow on such a bounded domain
obtained with the initial condition (3.34) develops a singularity at the
origin after a finite, albeit quite long, time.

####  3.3.3 Diagnostic quantities

We analyze numerical solutions of system (1.1) based on the following
diagnostic quantities: the $\dot{H}^{3}$ seminorm of the velocity field
$\bm{u}(t)$, the $L^{\infty}$ norm of the vorticity field $\bm{\omega}(t)$ and
the width $\delta(t)$ of the analyticity strip characterizing the velocity
field, all of which are functions of time $t$.

We approximate $\|\bm{u}(t)\|_{\dot{H}^{3}}$ and
$\|\bm{\omega}(t)\|_{L^{\infty}}$ as

$\displaystyle\|\bm{u}(t)\|_{\dot{H}^{3}}$
$\displaystyle\approx\sqrt{\sum_{j_{1}=0}^{N-1}\sum_{j_{2}=0}^{N-1}\sum_{j_{3}=0}^{N-1}|2\pi\bm{j}|^{6}\left|\hat{\bm{u}}_{\bm{j}}(t)\right|^{2}},$
(3.36a) $\displaystyle\|\bm{\omega}(t)\|_{L^{\infty}}$
$\displaystyle\approx\max_{0\leq
j_{1},j_{2},j_{3}<N}\left|\bm{\omega}\left(j_{1}/N,j_{2}/N,j_{3}/N,t\right)\right|.$
(3.36b)

The width of the analyticity strip $\delta(t)$ has been used to diagnose
singularity formation in solutions of various PDEs (Sulem et al., 1983),
including the 3D Euler equations (1.1) (Bustamante & Brachet, 2012). To
approximate it, we first define the energy spectrum of $\bm{u}(\bm{x},t)$ as

$e(k,t):=\frac{1}{2}\sum_{j\leq|\bm{j}|<j+1}|\hat{\bm{u}}_{\bm{j}}(t)|^{2},\qquad
k=2\pi j,\quad j\in\mathbb{N},$ (3.37)

such that the kinetic energy of the flow is

$E(t):=\frac{1}{2}\int_{\mathbb{T}^{3}}|\bm{u}(\bm{x},t)|^{2}\,d\bm{x}=\frac{1}{2}\sum_{\bm{j}\in\mathbb{Z}^{3}}|\hat{\bm{u}}_{\bm{j}}(t)|^{2}=\sum_{j}e(2\pi
j,t).$ (3.38)

When $\bm{u}(t)$ is real-analytic and admits an extension to the complex space
$\mathbb{C}^{3}$, it can be shown that the energy spectrum (3.37) has the
following asymptotic representation as $k\rightarrow\infty$ (Carrier et al.,
2005)

$e(k,t)=C(t)k^{-n(t)}e^{-2k\delta(t)},$ (3.39)

where $C(t)$ is a scaling constant and $n(t)$ is the order of the singularity
nearest to the real line. Clearly, the loss of analyticity and hence the
singularity formation are signalled by the vanishing of $\delta(t)$. Fixing
the time $t$ and computing the logarithm of (3.39), the parameters $C(t)$,
$n(t)$ and $\delta(t)$ can be determined by performing a least-squares fit to
minimize the error functional

$\chi(t):=\sum_{k=k_{<}}^{k_{>}}\left[\ln e(k,t)-\ln C(t)+n(t)\ln
k+2k\delta(t)\right]^{2},$ (3.40)

where $k_{<}\geq 2\pi$, $k_{>}<2\pi\lfloor N/3\rfloor$.

##  4 Numerical results

To distinguish between regular and potentially singular evolution, we solve
Problem 2.2 on “short” and “long” time intervals $[0,T]$ and perform a
resolution refinement, where optimal initial conditions
$\widetilde{\bm{\eta}}_{T}^{N}$ and the corresponding Euler flows
$\bm{u}^{N}\left(t;\widetilde{\bm{\eta}}_{T}^{N}\right)$ are computed with
increasing numerical resolutions $N^{3}$, cf. § 3.3.1. When $T$ is
sufficiently small, such that the Euler system (1.1) is guaranteed by Theorem
1.1 to be well-posed on $[0,T]$, we expect that

$\lim_{N\to\infty}\Phi_{T}\left(\widetilde{\bm{\eta}}_{T}^{N}\right)<\infty,$
(4.1)

i.e., the objective functional (2.7) must remain finite upon resolution
refinement.

On the other hand, when $T$ is large enough such that Theorem 1.1 does not
assert the existence of classical solutions on $[0,T]$ for regular initial
data $\bm{\eta}$, the divergence of the objective functional (2.7) upon
resolution refinement

$\lim_{N\to\infty}\Phi_{T}\left(\widetilde{\bm{\eta}}_{T}^{N}\right)=\infty$
(4.2)

would signal the possibility of a singularity formation at some $t\in[0,T]$.
As we shall see in § 4.1 and § 4.2 below, behaviors consistent with (4.1) and
(4.2) are in fact observed in solutions of Problem 2.2 obtained with $T=25$
and $T=75$, respectively, with the caveat that the limits in (4.1)–(4.2) are
extrapolated from four resolutions only
($N^{3}=128^{3},256^{3},512^{3},1024^{3}$).

###  4.1 Results for $T=25$

Since solutions of Problem 2.2 depend on the value of the parameter $\sigma$
defining the Gevrey space via (2.2) and (2.4), we begin by examining the
effect of this parameter on these solutions. We solve Problem 2.2 in different
Gevrey spaces $G^{\sigma}$ with $\sigma=10^{-1},10^{-2},\ldots,10^{-5}$ using
the resolution $128^{3}$ and two initial guesses $\bm{\eta}_{\text{TG}}$ and
$\bm{\eta}_{\text{rand}}$. The dependence of the value of the objective
functional $\Phi_{25}\left(\bm{\eta}^{(n)}\right)$ on the iteration index $n$
in (3.20) is shown for different $\sigma$ in figure 2 for the two initial
guesses. We observe that as long as $\sigma$ is sufficiently small
($\sigma<10^{-1}$), its value does not have an appreciable effect on the
maximum attained value $\widetilde{\Phi}_{25}$ of the objective functional.
However, this parameter does affect the rate of convergence and smaller values
of $\sigma$ typically result in faster convergence. This can be interpreted in
terms of an “effective dimension” of the discrete space used to approximate
the optimal initial condition $\widetilde{\bm{\eta}}_{25}^{N}$, which is
larger for smaller $\sigma$ since decreasing the value of this parameter
reduces the penalty on the Fourier components of
$\widetilde{\bm{\eta}}_{25}^{N}$ with large wavenumbers $|\bm{k}|$. Hereafter
we will use $\sigma=10^{-5}$ which results in both larger attained values of
$\widetilde{\Phi}_{25}$ and faster convergence of iterations in (3.20), cf.
Figure 2.

Figure 2: [Short time window, $T=25$] Dependence of the objective functional
$\Phi_{25}\left(\bm{\eta}^{(n)}\right)$ on the iteration index $n$ for
different values of $\sigma$ when (a) $\bm{\eta}_{{\text{TG}}}$ and (b)
$\bm{\eta}_{{\text{rand}}}$ are used as initial guesses in (3.20).

Since Problem 2.2 is non-convex, iterations (3.20) may lead to different local
maximizers depending on the choice of the initial guess $\bm{u}^{0}$. To
assess this possibility, we solve Problem 2.2 using four initial guesses
$\bm{\eta}_{\text{TG}}$, $\bm{\eta}_{\text{rand}}$, $\bm{\eta}_{\text{K}}$ and
$\bm{\eta}_{\text{H}}$ introduced in § 3.3.2 with the resolution $128^{3}$.
The results are presented in figures 3(a) and 3(b) where we show,
respectively, the value of the objective functional
$\Phi_{25}\left(\bm{\eta}^{(n)}\right)$ as a function of the iteration index
$n$ and the energy spectra (3.37) of the optimal initial conditions
$\widetilde{\bm{\eta}}_{25}$ obtained with different initial guesses. From the
coincidence of the terminal values $\widetilde{\Phi}_{25}$ in figure 3(a) and
of the energy spectra $e(k,0)$ in figure 3(b) we conclude that an essentially
the same optimal initial condition $\widetilde{\bm{\eta}}_{25}$ is found when
the initial guesses $\bm{\eta}_{\text{TG}}$, $\bm{\eta}_{\text{rand}}$ and
$\bm{\eta}_{\text{H}}$ are used. On the other hand, with the initial guess
$\bm{\eta}_{\text{K}}$, a much lower (approximately by a factor of 4) terminal
value $\widetilde{\Phi}_{25}$ of the objective functional is found, cf. figure
3(a), which corresponds to an optimal initial condition with less energy in
the high-wavenumber part of the spectrum in figure 3(b). The reason for this
is that the initial condition $\bm{\eta}_{\text{K}}$ is designed to produce a
large growth of different diagnostic quantities only after a long time,
approximately $t\approx 110059$ (this value of $t$ is obtained by rescaling
the time when Hou & Li (2006) stopped their computations to the normalization
$\|\bm{\eta}\|_{\dot{H}^{3}}=1$ we use in the present study). In the reminder
of this subsection we will use $\bm{\eta}_{\text{TG}}$ as the initial guess in
(3.20), which will allow us to compare our results with the findings of
earlier studies of potential singularities in Euler flows that used the
Taylor-Green vortex as the initial condition (Bustamante & Brachet, 2012).

Figure 3: [Short time window, $T=25$] (a) Dependence of the objective
functional $\Phi_{25}\left(\bm{\eta}^{(n)}\right)$ on the iteration index $n$
and (b) energy spectra (3.37) of the optimal initial conditions
$\widetilde{\bm{\eta}}_{25}$ corresponding to different initial guesses
$\bm{\eta}_{{\text{TG}}}$, $\bm{\eta}_{{\text{rand}}}$,
$\bm{\eta}_{{\text{K}}}$ and $\bm{\eta}_{{\text{H}}}$ used in iterations
(3.20). In panel (b) the solid line represents the Gaussian filter we use (Hou
& Li, 2007).

As is evident from figure 3(b), calculations performed with the resolution
$128^{3}$ are not fully resolved since in all cases the optimal initial
conditions $\widetilde{\bm{\eta}}_{25}$ have Fourier coefficients with
magnitudes larger than the machine precision for wavenumbers
$|\bm{k}|>k_{0}:=2\pi\lfloor N/3\rfloor$, which results in aliasing errors. To
address this issue, we solve Problem 2.2 with increasing resolutions
$N^{3}=256^{3},512^{3},1024^{3}$, which is done using the optimal initial
condition $\widetilde{\bm{\eta}}_{25}^{N}$ as the initial guess for the
iteration (3.20) performed with the resolution $(2N)^{3}$. The dependence of
the objective functional $\Phi_{25}^{N}\left(\bm{\eta}^{(n)}\right)$ on the
iteration index $n$ is shown in figure 4(a), where we see that larger values
of the objective functional are achieved each time the resolution is refined.
We also observe that after each resolution refinement, fewer iterations are
needed to achieve convergence; this is because only the Fourier coefficients
with high wavenumbers, which represent the fine structures in the flow, need
to be adjusted following a resolution refinement. As is evident from figure
4(b), the maximum values of the objective functional
$\widetilde{\Phi}_{25}^{N}$ converge to a well-defined finite limit as the
resolution increases, cf. (4.1).

The energy spectra of the optimal initial conditions
$\widetilde{\bm{\eta}}_{25}^{N}$ and of the corresponding terminal states
$\bm{u}^{N}\left(25;\widetilde{\bm{\eta}}_{25}^{N}\right)$ are shown in
figures 5(a) and 5(b) for different resolutions. These plots indicate that the
computations are already well resolved with $N^{3}=512^{3}$. The time
evolution of the width $\delta(t)$ of the analyticity strip and of the
corresponding order $n(t)$ of the singularity are shown for
$N^{3}=512^{3},1024^{3}$ in figures 6(a) and 6(b), respectively. In the first
figure we see that $\delta(t)$ reveals a mild decrease only near the end of
the time window $[0,25]$. These observations further confirm that the extreme
Euler flows which correspond to the initial data $\widetilde{\bm{\eta}}_{25}$
obtained as local maximizers of Problem 2.2 for $T=25$ remain regular on the
interval $[0,25]$.

Figure 4: [Short time window, $T=25$] Dependence of (a) the objective
functional $\Phi_{25}^{N}\left(\bm{\eta}^{(n)}\right)$ on the iteration index
$n$ for different resolutions $N^{3}$ and (b) of the corresponding maximum
attained values $\widetilde{\Phi}_{25}^{N}$ of the objective functional on
$N$.

Figure 5: [Short time window, $T=25$] The energy spectra of (a) the optimal
initial conditions $\widetilde{\bm{\eta}}_{25}^{N}$ and (b) of the
corresponding terminal states
$\bm{u}^{N}\left(25;\widetilde{\bm{\eta}}_{25}^{N}\right)$ obtained for
different resolutions $N^{3}$. The solid lines represent the Gaussian filters
we use (Hou & Li, 2007) whereas the dashed lines mark the threshold wavenumber
$k_{0}$ above which aliasing errors occur.

To close this subsection, in figures 7(a) and 7(b) we compare time evolution
of the norms $\|\bm{u}(t)\|_{\dot{H}^{3}}$ and
$\|\bm{\omega}(t)\|_{L^{\infty}}$ in Euler flows corresponding to different
initial conditions listed in § 3.3.2. We see that the growth of
$\|\bm{u}(t)\|_{\dot{H}^{3}}$ is by far the largest in the flow with the
optimal initial condition $\widetilde{\bm{\eta}}_{25}^{1024}$ and is only
modest when the other initial conditions are used. Likewise, the norm
$\|\bm{\omega}(t)\|_{L^{\infty}}$ exhibits some weak growth only in the flow
with the optimal initial condition $\widetilde{\bm{\eta}}_{25}^{1024}$ while
it actually decreases for $t\in[0,25]$ when the initial conditions
$\bm{\eta}_{{\text{TG}}}$, $\bm{\eta}_{{\text{rand}}}$ and
$\bm{\eta}_{\text{H}}$ are used.

Figure 6: [Short time window, $T=25$] Dependence of (a) the width $\delta(t)$
of the analyticity strip and (b) the corresponding order $n(t)$ of the
singularity, cf. (3.39), on time $t\in[0,25]$ in the extreme flows computed
with resolutions $N^{3}=512^{3},1024^{3}$.

Figure 7: [Short time window, $T=25$] Dependence of the solution norms (a)
$\|\bm{u}(t)\|_{\dot{H}^{3}}$ and (b) $\|\bm{\omega}(t)\|_{L^{\infty}}$ for
$t\in[0,25]$ in Euler flows corresponding to different initial conditions and
approximated using the resolution $N^{3}=1024^{3}$.

###  4.2 Results for $T=75$

We now move on to discuss solutions of Problem 2.2 on a longer time window
with $T=75$ where the behavior of the maximizing flows is qualitatively
different from their behavior on the short time window discussed in § 4.1. We
begin our discussion by analyzing the effect of the initial guess
$\bm{\eta}_{0}$ on the local maximizers of Problem 2.2 obtained with the
iterative algorithm (3.20). In figure 8(a) we show the dependence of the
objective functional $\Phi_{75}\left(\bm{\eta}^{(n)}\right)$ on the iteration
index $n$ and figure 8(b) shows the energy spectra of the corresponding
optimal initial conditions $\widetilde{\bm{\eta}}_{75}$. We remark that, as an
intermediate step, we first use the optimal initial condition
$\widetilde{\bm{\eta}}_{25}^{128}$ as the initial guess to solve Problem 2.2
with $T=50$, and then use its solution, $\widetilde{\bm{\eta}}_{50}^{128}$, as
the initial guess for $T=75$. We see that, with the exception of the initial
guess $\bm{\eta}_{\text{K}}$, computations performed with all other initial
guesses lead to the same local maximizer (up to rotations and translations).
As was the case when Problem 2.2 was solved with $T=25$, computations using
the initial guess $\bm{\eta}_{\text{K}}$ lead to a much smaller maximum
attained value of the objective functional since this initial condition is
designed to promote a significant growth of various regularity indicators,
such as $\|\bm{u}(t)\|_{\dot{H}^{3}}$ and $||\bm{\omega}(t)||_{L^{\infty}}$,
only on much longer time scales (Hou & Li, 2006).

Figure 8: [Long time window, $T=75$] (a) Dependence of the objective
functional $\Phi_{75}\left(\bm{\eta}^{(n)}\right)$ on the iteration index $n$
and (b) energy spectra (3.37) of the optimal initial conditions
$\widetilde{\bm{\eta}}_{75}$ corresponding to different initial guesses
$\bm{\eta}_{{\text{TG}}}$, $\bm{\eta}_{{\text{rand}}}$,
$\bm{\eta}_{{\text{K}}}$ and $\widetilde{\bm{\eta}}_{25}$ used in iterations
(3.20). In panel (b) the solid line represents the Gaussian filter we use (Hou
& Li, 2007).

As is evident from figure 8(b), computations carried out with the resolution
$128^{3}$ are under-resolved and now we consider the effect of refining the
resolution on solutions of Problem 2.2 with $T=75$. As we did in § 4.1, we
proceed by using the optimal initial condition
$\widetilde{\bm{\eta}}_{75}^{N}$ as the initial guess in iteration (3.20) when
solving Problem 2.2 with the resolution $(2N)^{3}$. The dependence of the
objective functional $\Phi_{75}^{N}\left(\bm{\eta}^{(n)}\right)$ on the
iteration index $n$ is presented for different resolutions in figure 9(a).
Figure 9(b) shows that, unlike in the case with $T=25$, cf. figure 4(b), the
maximum attained values of the objective functional
$\widetilde{\Phi}_{75}^{N}$ diverge as we refine the resolution with the
difference $\widetilde{\Phi}_{75}^{2N}-\widetilde{\Phi}_{75}^{N}$ increasing
with $N$. As explained in § 4, this indicates the possibility of a singularity
formation at some $t\in[0,75]$.

Figure 9: [Long time window, $T=75$] Dependence of (a) the objective
functional $\Phi_{75}^{N}\left(\bm{\eta}^{(n)}\right)$ on the iteration index
$n$ for different resolutions $N^{3}$ and (b) of the corresponding maximum
attained values $\widetilde{\Phi}_{75}^{N}$ of the objective functional on
$N$.

The energy spectra of the optimal initial conditions
$\widetilde{\bm{\eta}}_{75}^{N}$ and the corresponding terminal states
$\bm{u}^{N}\left(75;\widetilde{\bm{\eta}}_{75}^{N}\right)$ are shown in
figures 10(a) and 10(b), respectively. Comparing these two figures, we observe
that even though the initial conditions $\widetilde{\bm{\eta}}_{75}^{N}$ are
well-resolved for $N=512,1024$, the time evolutions computed with all
resolutions become under-resolved at later times. In particular, for $N=1024$
this happens at $t=51.25$ and is detected by the appearance of aliasing errors
when the amplitudes of Fourier coefficients with wavenumbers
$|\bm{k}|>2\pi\lfloor N/3\rfloor$ become larger than the machine precision.
For lower resolutions, computations become under-resolved at earlier times.

Figure 10: [Long time window, $T=75$] The energy spectra of (a) the optimal
initial conditions $\widetilde{\bm{\eta}}_{75}^{N}$ and (b) of the
corresponding terminal states
$\bm{u}^{N}\left(75;\widetilde{\bm{\eta}}_{75}^{N}\right)$ obtained for
different resolutions $N^{3}$. The solid lines represent the Gaussian filters
we use (Hou & Li, 2007) whereas the dashed lines mark the threshold wavenumber
$k_{0}$ above which aliasing errors occur.

The time evolution of the width $\delta(t)$ of the analyticity strip and of
the corresponding order $n(t)$ of the singularity in optimal flows computed
with $N=512,1024$ is shown in figures 11(a) and 11(b), respectively. These
quantities are shown only up to times when they can still be determined in a
reliable manner, i.e., as long as $e(2\pi\lfloor N/3\rfloor,t)\leq 10^{-20}$.
In figure 11(a) we see that, following an initial period of growth,
$\delta(t)$ starts to decrease around the middle of the time window and,
interestingly, this decrease becomes faster than exponential approximately at
the time when the computations become under-resolved.

Figure 11: [Long time window, $T=75$] Dependence of (a) the width $\delta(t)$
of the analyticity strip and (b) the corresponding order $n(t)$ of the
singularity, cf. (3.39), on time $t\in[0,75]$ in the extreme flows computed
with resolutions $N^{3}=512^{3},1024^{3}$.

The time evolution of $\|\bm{u}(t)\|_{\dot{H}^{3}}$ and
$\|\bm{\omega}(t)\|_{L^{\infty}}$ in the flows corresponding to different
initial conditions is shown in figures 12(a) and 12(b), respectively. As was
the case in § 4.1, a significant growth of these norms is evident only in the
flows corresponding to the optimal initial condition
$\widetilde{\bm{\eta}}_{75}^{1024}$, whereas $\|\bm{\omega}(t)\|_{L^{\infty}}$
actually decreases on $[0,75]$ for solutions with the initial conditions
$\bm{\eta}_{{\text{TG}}}$, $\bm{\eta}_{{\text{rand}}}$ and
$\bm{\eta}_{{\text{H}}}$.

Figure 12: [Long time window, $T=75$] Dependence of the solution norms (a)
$\|\bm{u}(t)\|_{\dot{H}^{3}}$ and (b) $\|\bm{\omega}(t)\|_{L^{\infty}}$ for
$t\in[0,75]$ in Euler flows corresponding to different initial conditions and
approximated using the resolution $N^{3}=1024^{3}$.

In order to shed light on whether the growth of the norm
$\|\bm{u}(t)\|_{\dot{H}^{3}}$ observed in figure 12(a) corresponding to the
optimal initial conditions $\widetilde{\bm{\eta}}_{75}^{N}$ may indicate a
singularity formation at some time $t\in[0,75]$, we further analyze the growth
rate of $\|\bm{u}(t)\|_{\dot{H}^{3}}$. We assume that the evolution of this
norm is described by the relation

$\frac{d\|\bm{u}(t)\|_{\dot{H}^{3}}}{dt}=C(t)\|\bm{u}(t)\|_{\dot{H}^{3}}^{\alpha(t)},$
(4.3)

which is motivated by the structure of rigorous a priori estimates for the
rate of growth of Sobolev norms of solutions to the Navier-Stokes system
(Doering & Gibbon, 1995; Lu & Doering, 2008; Ayala & Protas, 2017). We note
that if the exponent $\alpha(t)$ in (4.3) remains larger than 1 over a
sufficiently long time with the prefactor $C(t)$ bounded away from zero, then
there will be a finite-time blow-up of the norm $\|\bm{u}(t)\|_{\dot{H}^{3}}$,
and thus formation of a singularity in the corresponding Euler flow, cf.
Theorem 1.1. To obtain further insights about these quantities, in figure
13(a) we plot $(d/dt)\|\bm{u}(t)\|_{\dot{H}^{3}}$ versus
$\|\bm{u}(t)\|_{\dot{H}^{3}}$ with $t\in[0,75]$ using log-log scaling for the
optimal flows computed with different resolutions, such that the exponent
$\alpha(t)$ can be inferred from the slope of the tangent to the curves at
$\|\bm{u}(t)\|_{\dot{H}^{3}}$. The evolution of the exponent $\alpha(t)$
determined by a local fitting procedure applied to ansatz (4.3) with time
$t\in[0,75]$ is shown for $N=1024$ in figure 13(b), where a decreasing trend
is evident and we have $\alpha(t)>1$ for $t\in[0,53.495]$ while the
computation becomes under-resolved at $t\approx 51.25$. As regards the
prefactor $C(t)$ in (4.3), it reveals a slow growth with time $t$ (and with
$\|\bm{u}(t)\|_{\dot{H}^{3}}$) which is well approximated by the expression
$C(t)=0.0568\left(\ln\|\bm{u}(t)\|_{\dot{H}^{3}}\right)^{0.5742}$ with
parameters determined via a least-squares fit. We thus conclude that the time
evolution of the norm $\|\bm{u}(t)\|_{\dot{H}^{3}}$ in the flow with the
optimal initial condition $\widetilde{\bm{\eta}}_{75}^{1024}$ remains
consistent with the singularity formation as long as the computation remains
well-resolved.

Figure 13: [Long time window, $T=75$] (a) Dependence of
$(d/dt)\|\bm{u}(t)\|_{\dot{H}^{3}}$ on $\|\bm{u}(t)\|_{\dot{H}^{3}}$ in Euler
flows with the optimal initial conditions $\widetilde{\bm{\eta}}_{75}^{N}$ for
$t\in[0,75]$ and different resolutions $N^{3}$ and (b) the corresponding
exponent $\alpha(t)$ in (4.3) for $t\in[0,75]$ with $N=1024$. The straight
lines in both panels represent $\alpha=1$.

Finally, we analyze the physical-space structure of the extreme flow,
beginning with the vorticity field of the optimal initial condition
$\widetilde{\bm{\omega}}_{75}^{1024}=[\omega_{1},\omega_{2},\omega_{3}]=\bm{\nabla}\times\widetilde{\bm{\eta}}_{75}^{1024}$.
In figure 14, we show the three components of
$\widetilde{\bm{\omega}}_{75}^{1024}$ and observe that the optimal initial
condition has the form of three perpendicular pairs of distorted anti-parallel
vortex tubes. If we define the component-wise enstrophy as
$\mathcal{E}_{i}:=\frac{1}{2}\int_{\mathbb{T}^{3}}\omega_{i}^{2}(\bm{x})\;d\bm{x}$,
$i=1,2,3$, then we have $\mathcal{E}_{1}=6.34\times
10^{-5}\approx\mathcal{E}_{2}=6.35\times 10^{-5}<\mathcal{E}_{3}=8.03\times
10^{-5}$, which is a signature of the symmetry of the flow with respect to the
plane $x_{1}=x_{2}$. In fact, this initial condition possesses a similar
physical-space structure as the one found in Kang et al. (2020), where the
authors searched for singularities in Navier-Stokes flows by maximizing the
total enstrophy at a prescribed final time. While, unlike here, no evidence
for singularity formation was revealed in that earlier study.

Next, in figure 15(a), we analyze the optimal initial condition
$\widetilde{\bm{\eta}}_{75}^{1024}$ in more detail with the corresponding
terminal state
$\bm{u}^{1024}\left(75;\widetilde{\bm{\eta}}_{75}^{1024}\right)$ shown in
figure 15(b). In both figures we show the isosurfaces of $|\bm{\omega}|$ and
$\log_{10}\left(\left||D|^{3}\bm{u}\right|\right)$ together with selected
streamlines passing through a small neighborhood of the origin (see also the
animation available at https://youtu.be/G_nfJNq6W1A). We observe that
throughout its evolution the optimal flow has the form of two jets colliding
head-on with the vorticity concentrating into two strongly flattened vortex
rings as time goes on. The region with large values of
$\log_{10}\left(\left||D|^{3}\bm{u}\right|\right)$, which is the quantity that
matters in our objective functional (2.7), evolves into a flat disc located in
the middle of the two rings. The vorticity field at the final time $t=75$ has
three symmetry planes: $x_{1}=x_{2}$, $x_{1}=x_{3}$ and $x_{2}=x_{3}$, in
addition to discrete rotation symmetries with respect to the body diagonal
passing through the center of the two rings. Without loss of generality, we
focus our discussion on the symmetry plane $x_{1}=x_{2}$, and in figures 16(a)
and 16(b) we visualize the vorticity component
$\omega^{\perp}:=\bm{\omega}\cdot\bm{n}$ normal to that plane
($\bm{n}=\left[-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\right]^{T}$ is the
unit vector normal to the symmetry plane), at $t=0$ and $t=75$ (see also the
animation available at https://youtu.be/6BcQ4LvcAE4). The streamline pattern
in Figure 15 (b) indicates that as $t\rightarrow 75$ the flow in the two jets
colliding near the origin sharply transitions towards a radial outflow through
the gap between the two vortex rings. At the same time, this gap becomes very
narrow, which results in a sharp transition between the regions of the
symmetry plane characterized by opposite signs of the normal vorticity
$\omega^{\perp}$, cf. figure 16(b). This appears to be the mechanism
responsible for the possible singularity formation in the extreme flow
considered here.

Figure 14: [Long time window, $T=75$] Isosurfaces of the components (a)
$\omega_{1}$, (b) $\omega_{2}$ and (c) $\omega_{3}$ of the vorticity field
corresponding to the optimal initial condition
$\widetilde{\bm{\eta}}_{75}^{1024}$.

Figure 15: [Long time window, $T=75$] Isosurfaces of $|\bm{\omega}|$ and
$\log_{10}(||D|^{3}\bm{u}|)$ in (a) the optimal initial condition
$\widetilde{\bm{\eta}}_{75}^{1024}$ and (b) the corresponding terminal state
$\bm{u}^{1024}\left(75;\widetilde{\bm{\eta}}_{75}^{1024}\right)$ together with
selected streamlines. The same color ranges are used in both panels. An
animated version of these figures (without streamlines) showing the time
evolution for $t\in[0,75]$ is available https://youtu.be/G_nfJNq6W1A.

Figure 16: [Long time window, $T=75$] The vorticity component $\omega^{\perp}$
normal to the symmetry plane $x_{1}=x_{2}$ in (a) the optimal initial
condition $\widetilde{\bm{\eta}}_{75}^{1024}$ and (b) the corresponding
terminal state
$\bm{u}^{1024}\left(75;\widetilde{\bm{\eta}}_{75}^{1024}\right)$. The same
color range is used in both panels. An animated version of these figures
showing the time evolution for $t\in[0,75]$ is available at
https://youtu.be/6BcQ4LvcAE4.

##  5 Discussion and conclusions

In this paper, we study the problem concerning the possibility of spontaneous
formation of singularities in solutions of the 3D Euler equations (1.1) on a
periodic domain. Based on the local well-posedness results by Kato (1972)
stated in Theorem 1.1, we formulate a PDE-constrained optimization problem to
search for initial conditions $\bm{\eta}$ with unit $\dot{H}^{3}$ seminorm
such that the corresponding optimal solution achieves a maximum $\dot{H}^{3}$
seminorm at a prescribed time $T$. Since we focus on smooth (real-analytic)
initial data, the optimization problem is formulated in a suitable Gevrey
space (2.4). It is then solved in the “optimize-then-discretize” setting using
a state-of-the-art Riemannian conjugate gradient method (3.20) where the
search direction at every iteration depends on the gradient of the objective
functional at the current iteration as well as on the previous search
direction; the former is obtained by solving the adjoint system (3.8) backward
in time. The required regularity of the gradient is ensured by the use of the
Riesz representation theorem, cf. (3.4). An analogous approach has been
successfully used to solve PDE-constrained optimization problems formulated to
elucidate extreme behaviors in 1D Burgers (Ayala & Protas, 2011) and 3D
Navier-Stokes flows (Kang et al., 2020; Kang & Protas, 2021).

Problem 2.2 is non-convex and we have found evidence for the presence of
nonunique local maximizers corresponding to different initial guesses listed
in § 3.3.2. However, iterations performed with most initial guesses were found
to converge to the same (up to rotation and translation) local maximizer, with
the exception of the initial guess $\bm{\eta}_{\text{K}}$, which was designed
to produce significant growth of various regularity indicators at times much
longer than the intervals $[0,T]$ considered in this study.

We adopt an indirect approach to determine whether a singularity may occur
within the time interval $[0,T]$ based on sequential refinements of the
numerical resolution. Specifically, we solve Problem 2.2 with $T=25$ and
$T=75$ using spatial resolutions increasing from $128^{3}$ to $1024^{3}$. For
$T=25$, we observe that the maximum attained values of
$\left|\left|\bm{u}^{N}\left(75;\widetilde{\bm{\eta}}_{75}^{N}\right)\right|\right|_{\dot{H}^{3}}$
converge to a finite limit for the optimal solutions obtained with increasing
resolutions, which indicates that the Euler system (1.1) is well-posed on
$[0,25]$. However, for $T=75$, the maximum values of
$\left|\left|\bm{u}^{N}\left(75;\widetilde{\bm{\eta}}_{75}^{N}\right)\right|\right|_{\dot{H}^{3}}$
attained in the optimal solutions diverge upon resolution refinement, which
suggests a possible formation of singularity for some $t\in[0,75]$. As shown
in figure 13, the growth of the norm $\|\bm{u}(t)\|_{\dot{H}^{3}}$ in the
flows corresponding to the optimal initial condition obtained for $T=75$ using
the highest numerical resolution proceeds at a rate consistent with a finite-
time blow-up for as long as the solution remains well-resolved, although this
rate slows down with time. With the limited numerical resolution we could use,
it is impossible to conclude whether this depletion of the growth rate would
eventually prevent a singularity from appearing in a finite time. We add that
the largest resolution we used ($1024^{3}$) is smaller than the largest
resolutions available to-date in simulations of the Navier-Stokes and the
Euler systems. However, we emphasize that in the present study we solve a
family of optimization problems which are much more costly to solve than
simply computing the evolution of the solution on the same time interval. In
fact, solution of Problem 2.2 typically requires $\mathcal{O}(100)$ solutions
of the Euler system (1.1) and $\mathcal{O}(10)$ solutions of the adjoint
system (3.8).

We add that under-resolved computations of Euler flows lead to
“thermalization” of solutions, a process where after a sufficiently long time
the kinetic energy is equidistributed among the finite number of Fourier modes
used in the computation. Thermalization is related to the appearance of the
so-called “tygers” and inevitably occurs when using numerical methods based on
Galerkin truncation to compute the time evolution of inviscid systems, such as
1D Burgers equation (Rampf et al., 2022) and 3D incompressible Euler equations
(Murugan & Ray, 2022).

Figures 7 and 12 show that, as compared to flows obtained with different
initial guesses, the optimal flows demonstrate a more significant growth in
both $\|\bm{u}(t;\widetilde{\bm{\eta}})\|_{\dot{H}^{3}}$ and
$\|\bm{\omega}(t;\widetilde{\bm{\eta}})\|_{L^{\infty}}$ throughout the
optimization time windows $[0,T]$ with $T=25$ and $T=75$. We further notice
that the growth of the norm $\|\bm{\omega}(t)\|_{L^{\infty}}$ in the extreme
flows obtained by solving Problem 2.2 with $T=75$ is weaker than the growth of
the norm $\|\bm{u}(t)\|_{\dot{H}^{3}}$ in the same flow, cf. figures 12(a) and
12(b). However, estimate (1.4) shows that an unbounded growth of the latter
norm implies such a growth for the former norm as well and the doubly-
exponential structure of the upper bound in this estimate explains why in a
blow-up scenario the growth of $\|\bm{\omega}(t)\|_{L^{\infty}}$ may be weaker
than that of $\|\bm{u}(t)\|_{\dot{H}^{3}}$. As regards the behavior of the
width $\delta(t)$ of the analyticity strip in these flows, we note that the
decay rate of this quantity accelerates and becomes faster than exponential
around the time when computations become under-resolved.

In regard to the structure of the extreme flows in the physical space, it is
interesting to see in figure 15(b) that some key features of these flows are
close to being axisymmetric, to the extent that this may be possible for a
flow evolving on a cubic domain. More specifically, the region where the
largest values of $\left||D|^{3}\bm{u}\right|$ are concentrated has the form
of a circular disc embedded between the two colliding vortex rings, which
becomes more irregular further away from the axis of symmetry. We emphasize
that the extreme flows we found here are fairly robust in the sense that they
have been found consistently using several different initial guesses in
(3.20), albeit not with all initial guesses we used.

As concerns future work, the properties of the optimal flows found in the
present study suggest looking for extreme behaviors among initial conditions
constrained by some symmetries, such as axial and/or reflection symmetries, as
was pursued by Luo & Hou (2014a, b); Hou (2022). In addition, given that
spectral methods have a limited ability to resolve localized small-scale
features of potentially singular flows, further progress will likely require
the use of adaptive discretization techniques.

## Acknowledgements

The authors wish to express thanks to Marc Brachet, Miguel Bustamante, John
Gibbon, Robert Kerr, Adam Larios and Thomas Y. Hou for enlightening and
enjoyable discussions.The authors also acknowledge the support through an
NSERC (Canada) Discovery Grant. Computational resources were provided by
Digital Research Alliance of Canada under its Resource Allocation Competition.

## References

  * Absil et al. (2008) Absil, P.-A., Mahony, R. & Sepulchre, R. 2008 Optimization Algorithms on Matrix Manifolds. Princeton University Press.
  * Adams & Fournier (2005) Adams, R. A. & Fournier, J. F. 2005 Sobolev Spaces. Elsevier.
  * Ayala (2014) Ayala, D. 2014 Extreme Vortex States and Singularity Formation in Incompressible Flows. PhD thesis, McMaster University, available at http://hdl.handle.net/11375/15453.
  * Ayala & Protas (2011) Ayala, D. & Protas, B. 2011 On maximum enstrophy growth in a hydrodynamic system. Physica D 240, 1553–1563.
  * Ayala & Protas (2014) Ayala, D. & Protas, B. 2014 Maximum palinstrophy growth in 2D incompressible flows. Journal of Fluid Mechanics 742, 340–367.
  * Ayala & Protas (2017) Ayala, D. & Protas, B. 2017 Extreme vortex states and the growth of enstrophy in 3D incompressible flows. Journal of Fluid Mechanics 818, 772–806.
  * Bardos & Benachour (1977) Bardos, Claude & Benachour, S 1977 Domaine d’analycité des solutions de l’équation d’Euler dans un ouvert de $R^{n}$. Annali della Scuola Normale Superiore di Pisa-Classe di Scienze 4 (4), 647–687.
  * Beale et al. (1984) Beale, J. T., Kato, T. & Majda, A. 1984 Remarks on the breakdown of smooth solutions for the $3$-D Euler equations. Comm. Math. Phys. 94 (1), 61–66.
  * Brachet (1991) Brachet, M. E. 1991 Direct simulation of three-dimensional turbulence in the Taylor-Green vortex. Fluid Dynamics Research 8, 1–8.
  * Brachet et al. (1983) Brachet, M. E., Meiron, D. I., Orszag, S. A., Nickel, B. G., Morf, R. H. & Frisch, U. 1983 Small-scale structure of the Taylor-Green vortex. Journal of Fluid Mechanics 130, 411–452.
  * Bustamante & Brachet (2012) Bustamante, M. D. & Brachet, M. 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.
  * Bustamante & Kerr (2008) Bustamante, M. D. & Kerr, R. M. 2008 3D Euler about a 2D symmetry plane. Physica D 237, 1912–1920.
  * Campolina & Mailybaev (2018) Campolina, Ciro S. & Mailybaev, Alexei A. 2018 Chaotic Blowup in the 3D Incompressible Euler Equations on a Logarithmic Lattice. Phys. Rev. Lett. 121, 064501.
  * Carrier et al. (2005) Carrier, George F., Krook, Max & Pearson, Carl E. 2005 Functions of a Complex Variable: Theory and Technique. Society for Industrial and Applied Mathematics, arXiv: https://epubs.siam.org/doi/pdf/10.1137/1.9780898719116.
  * Chae (2001) Chae, Dongho 2001 On the Well-Posedness of the Euler Equations in the Besov and the Triebel-Lizorkin Spaces . Tosio Kato’s Method and Principle for Evolution Equations in Mathematical Physics pp. 42–57.
  * Chen & Hou (2022) Chen, Jiajie & Hou, Thomas Y 2022 Stable nearly self-similar blowup of the 2D Boussinesq and 3D Euler euqations with smooth data. arXiv preprint arXiv:2210.07191 .
  * Cichowlas & Brachet (2005) Cichowlas, C. & Brachet, M. E. 2005 Evolution of complex singularities in Kida-Pelz and Taylor-Green inviscid flows. Fluid Dynamics Research 36, 239–248.
  * Constantin et al. (1994) Constantin, Peter, E, Weinan & Titi, Edriss S. 1994 Onsager’s conjecture on the energy conservation for solutions of Euler’s equation. Communications in Mathematical Physics 165 (1), 207 – 209.
  * Constantin et al. (1996) Constantin, Peter, Fefferman, Charles & Majda, Andrew J 1996 Geometric constraints on potentially singular solutions for the 3-D Euler equations. Communications in Partial Differential Equations 21 (3-4).
  * Danaila & Protas (2017) Danaila, Ionut & Protas, Bartosz 2017 Computation of ground states of the gross–pitaevskii functional via riemannian optimization. SIAM Journal on Scientific Computing 39 (6), B1102–B1129, arXiv: https://doi.org/10.1137/17M1121974.
  * De Lellis & Székelyhidi Jr (2019) De Lellis, Camillo & Székelyhidi Jr, László 2019 On turbulence and geometry: from nash to onsager. Notices Amer. Math. Soc. 66 (5), 677–685, cvgmt preprint.
  * Deng et al. (2005) Deng, Jian, Hou, Thomas Y & Yu, Xinwei 2005 Geometric properties and non-blowup of 3D incompressible Euler flow. Communications in Partial Difference Equations 30 (1-2), 225–243.
  * Doering (2009) Doering, C. R. 2009 The 3D Navier-Stokes problem. Annual Review of Fluid Mechanics 41, 109–128.
  * Doering & Gibbon (1995) Doering, C. R. & Gibbon, J. D. 1995 Applied Analysis of the Navier-Stokes Equations. Cambridge University Press.
  * Elgindi (2021) Elgindi, Tarek M 2021 Finite-time singularity formation for $C^{1,\alpha}$ solutions to the incompressible Euler equations on $\mathbb{R}^{3}$ . Annals of Mathematics 194 (3), 647–727.
  * Elgindi & Jeong (2019) Elgindi, Tarek M. & Jeong, In-Jee 2019 Finite-Time Singularity Formation for Strong Solutions to the Axi-symmetric 3D Euler Equations. Annals of PDE 5 (2), 16.
  * Fefferman (2000) Fefferman, C. L. 2000 Existence and smoothness of the Navier-Stokes equation. available at http://www.claymath.org/sites/default/files/navierstokes.pdf, Clay Millennium Prize Problem Description.
  * Fehn et al. (2022) Fehn, Niklas, Kronbichler, Martin, Munch, Peter & Wall, Wolfgang A 2022 Numerical evidence of anomalous energy dissipation in incompressible Euler flows: towards grid-converged results for the inviscid Taylor–Green problem. Journal of Fluid Mechanics 932.
  * Frigo & Johnson (2003) Frigo, Matteo & Johnson, Steven G. 2003 FFTW User’s Manual. Massachusetts Institute of Technology.
  * Gibbon (2013) Gibbon, J.D. 2013 Dynamics of Scaled Norms of Vorticity for the Three-dimensional Navier-Stokes and Euler Equations. Procedia IUTAM 7, 39–48, iUTAM Symposium on Topological Fluid Dynamics: Theory and Applications.
  * Gibbon et al. (2008) Gibbon, J. D., Bustamante, M. & Kerr, R. M. 2008 The three–dimensional Euler equations: singular or non–singular? Nonlinearity 21, 123–129.
  * Grafke et al. (2008) Grafke, T., Homann, H., Dreher, J. & Grauer, R. 2008 Numerical simulations of possible finite-time singularities in the incompressible Euler equations: comparison of numerical methods. Physica D 237, 1932–1936.
  * Gunzburger (2003) Gunzburger, M. D. 2003 Perspectives in Flow Control and Optimization. SIAM.
  * Guo et al. (2022) Guo, Hui, Liang, Xueting & Yang, Yang 2022 Provable convergence of blow-up time of numerical approximations for a class of convection-diffusion equations. Journal of Computational Physics 466, 111421.
  * Hou (2009) Hou, T. Y. 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.
  * Hou (2022) Hou, Thomas Y 2022 Potential singularity of the 3D Euler equations in the interior domain. Foundations of Computational Mathematics pp. 1–47.
  * Hou & Li (2006) Hou, Thomas Y & Li, Ruo 2006 Dynamic depletion of vortex stretching and non-blowup of the 3-D incompressible Euler equations. Journal of Nonlinear Science 16 (6), 639–664.
  * Hou & Li (2007) Hou, T. Y. & Li, R. 2007 Computing nearly singular solutions using pseudo-spectral methods. Journal of Computational Physics 226, 379–397.
  * Kang & Protas (2021) Kang, D. & Protas, B. 2021 Searching for Singularities in Navier-Stokes Flows Based on the Ladyzhenskaya-Prodi-Serrin Conditions. ArXiv:2110.06130.
  * Kang et al. (2020) Kang, Di, Yun, Dongfang & Protas, Bartosz 2020 Maximum amplification of enstrophy in three-dimensional Navier-Stokes flows. Journal of Fluid Mechanics 893, A22.
  * Kato (1972) Kato, Tosio 1972 Nonstationary flows of viscous and ideal fluids in $\mathbb{R}^{3}$. Journal of functional Analysis 9 (3), 296–305.
  * Kerr (1993) Kerr, R. M. 1993 Evidence for a singularity of the three-dimensional, incompressible Euler equations. Phys. Fluids A 5, 1725–1746.
  * Kerr (2013) Kerr, Robert M. 2013 Swirling, turbulent vortex rings formed from a chain reaction of reconnection events. Physics of Fluids 25, 065101.
  * Kozono & Taniuchi (2000) Kozono, Hideo & Taniuchi, Yasushi 2000 Limiting Case of the Sobolev Inequality in BMO,with Application to the Euler Equations. Communications in Mathematical Physics 214 (1), 191–200.
  * Kukavica & Vicol (2009) Kukavica, Igor & Vicol, Vlad 2009 On the radius of analyticity of solutions to the three-dimensional Euler equations. Proceedings of the American Mathematical Society 137 (2), 669–677.
  * Lichtenstein (1925) Lichtenstein, Leon 1925 Über einige existenzprobleme der hydrodynamik homogener, unzusammendrückbarer, reibungsloser flüssigkeiten und die helmholtzschen wirbelsätze. Mathematische Zeitschrift 23 (1), 89–154.
  * Lu & Doering (2008) Lu, L. & Doering, C. R. 2008 Limits on enstrophy growth for solutions of the three-dimensional Navier–Stokes equations. Indiana University Mathematics Journal 57, 2693–2727.
  * Luenberger (1969) Luenberger, D. 1969 Optimization by Vector Space Methods. John Wiley and Sons.
  * Luo & Hou (2014a) Luo, G. & Hou, T. Y. 2014a Potentially Singular Solutions of the 3D Axisymmetric Euler Equations. Proceedings of the National Academy of Sciences 111 (36), 12968–12973.
  * Luo & Hou (2014b) Luo, G. & Hou, T. Y. 2014b Toward the Finite-Time Blowup of the 3D Incompressible Euler Equations: a Numerical Investigation. SIAM: Multiscale Modeling and Simulation 12 (4), 1722–1776.
  * Majda & Bertozzi (2002) Majda, A. J. & Bertozzi, A. L. 2002 Vorticity and Incompressible Flow. Cambridge University Press.
  * Matsumoto et al. (2008) Matsumoto, T., Bec, J. & Frisch, U. 2008 Complex-space singularities of 2D Euler flow in lagrangian coordinates. Physica D 237, 1951–1955.
  * Moffatt & Kimura (2019a) Moffatt, H. K. & Kimura, Yoshifumi 2019a Towards a finite-time singularity of the Navier-Stokes equations Part 1. Derivation and analysis of dynamical system. Journal of Fluid Mechanics 861, 930–967.
  * Moffatt & Kimura (2019b) Moffatt, H. K. & Kimura, Yoshifumi 2019b Towards a finite-time singularity of the Navier-Stokes equations. part 2. vortex reconnection and singularity evasion. Journal of Fluid Mechanics 870, R1.
  * Murugan & Ray (2022) Murugan, Sugan Durai & Ray, Samriddhi Sankar 2022 On the thermalization of the three-dimensional, incompressible, Galerkin-truncated Euler equation. arXiv preprint arXiv:2209.05046 .
  * Nocedal & Wright (1999) Nocedal, J. & Wright, S. J. 1999 Numerical Optimization. Springer.
  * Ohkitani (2008) Ohkitani, K. 2008 A miscellany of basic issues on incompressible fluid equations. Nonlinearity 21, 255–271.
  * Ohkitani & Constantin (2008) Ohkitani, K. & Constantin, P. 2008 Numerical study of the Eulerian–Lagrangian analysis of the Navier-Stokes turbulence. Phys. Fluids 20, 1–11.
  * Orlandi et al. (2014) Orlandi, P., Pirozzoli, S., Bernardini, M. & Carnevale, G. F. 2014 A minimal flow unit for the study of turbulence with passive scalars. Journal of Turbulence 15, 731–751.
  * Orlandi et al. (2012) Orlandi, P., Pirozzoli, S. & Carnevale, G. F. 2012 Vortex events in Euler and Navier-Stokes simulations with smooth initial conditions. Journal of Fluid Mechanics 690, 288–320.
  * Pelz (2001) Pelz, R. B. 2001 Symmetry and the hydrodynamic blow-up problem. Journal of Fluid Mechanics 444, 299–320.
  * Press et al. (1986) Press, W. H., Flannery, B. P., Teukolsky, S. A. & Vetterling, W. T. 1986 Numerical Recipes. Cambridge University Press.
  * Protas (2022) Protas, Bartosz 2022 Systematic search for extreme and singular behaviour in some fundamental models of fluid mechanics. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 380 (2225), 20210035, arXiv: https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2021.0035.
  * Pumir & Siggia (1990) Pumir, A. & Siggia, E. 1990 Collapsing solutions to the 3D Euler equations. Phys. Fluids A 2, 220–241.
  * Rampf et al. (2022) Rampf, Cornelius, Frisch, Uriel & Hahn, Oliver 2022 Eye of the tyger: early-time resonances and singularities in the inviscid burgers equation. arXiv preprint arXiv:2207.12416 .
  * Robinson (2020) Robinson, James C 2020 The Navier-Stokes regularity problem. Phil. Trans. R. Soc. A 378, 20190526.
  * Sato (2021) Sato, Hiroyuki 2021 Riemannian conjugate gradient methods: General framework and specific algorithms with convergence analyses, arXiv: 2112.02572.
  * Siegel & Caflisch (2009) Siegel, M. & Caflisch, R. E. 2009 Calculation of complex singular solutions to the 3D incompressible Euler equations. Physica D 238, 2368–2379.
  * Sulem et al. (1983) Sulem, C., Sulem, P. L. & Frisch, H. 1983 Tracing complex singularities with spectral methods. Journal of Computational Physics 50, 138–161.
  * Taylor & Green (1937) Taylor, G. I. & Green, A. E. 1937 Mechanism of the production of small eddies from large ones. Proceedings of the Royal Society of London A 158, 499–521.
  * Yin et al. (2021) Yin, Xi-Yuan, Schneider, Kai & Nave, Jean-Christophe 2021 A Characteristic Mapping Method for the three-dimensional incompressible Euler equations. arXiv preprint arXiv:2107.03504 .
  * Yun & Protas (2018) Yun, Dongfang & Protas, Bartosz 2018 Maximum Rate of Growth of Enstrophy in Solutions of the Fractional Burgers Equation. Journal of Nonlinear Science 28 (1), 395–422.