1. Introduction
Streamsurfaces of a three-dimensional flow are two-dimensional surfaces composed of streamlines. Even a small but judiciously chosen set of such surfaces can give an effective characterization of the global topology of the velocity field. In steady flows, streamsurfaces are also invariant manifolds for the particle motion and hence frame the Lagrangian particle dynamics. For these reasons, streamsurfaces should be, in principle, the simplest tool for illustrating instantaneous features of a velocity field.
Analogously, to identify and visualize vortical features of a velocity field, Yang & Pullin (Reference Yang and Pullin2010) generalized the notion of vortex tubes and sheets (Batchelor Reference Batchelor2000) by defining the vortex-surface field (VSF), i.e. a smooth scalar field whose isosurfaces act as two-dimensional (2-D) invariant manifolds of the vorticity field. Initially developed for symmetric, inviscid flows (Yang & Pullin Reference Yang and Pullin2010), this work has been extended to capture approximations to the Lagrangian (material) evolution of VSFs in analytic viscous flows (Yang & Pullin Reference Yang and Pullin2011), shear flows (Xiong & Yang Reference Xiong and Yang2017), compressible flows (Peng & Yang Reference Peng and Yang2018), transitional wall flows (Zhao, Yang & Chen Reference Zhao, Yang and Chen2016a,Reference Zhao, Yang and Chenb) and in homogeneous isotropic turbulence (Xiong & Yang Reference Xiong and Yang2019). A different approach, based on the spherical Clebsch maps, was used by Chern et al. (Reference Chern, Knöppel, Pinkall and Schröder2017) to visualize vortex lines and surfaces in computer graphics. However, neither VSFs nor streamsurfaces are objective (frame-indifferent; Haller Reference Haller2005), which renders their experimental detection ambiguous.
In contrast, recent results on the transport of dynamically active vector fields (such as the momentum and vorticity fields) have yielded objectively defined barrier vector fields whose distinguished invariant surfaces turn out to be frame-indifferent material barriers to active transport (Haller et al. Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020). Specifically, the barrier vector fields are the velocity and vorticity Laplacian in the instantaneous (Eulerian) limit and the time-averaged pullbacks of these Laplacian fields in the Lagrangian case. The invariant manifolds of these barrier vector fields have been shown to highlight vortical features of the velocity field with increased accuracy in several 2-D and three-dimensional (3-D) flows (Haller et al. Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020; Aksamit & Haller Reference Aksamit and Haller2022).
Irrespective of the underlying vector field, no general methodology is available for its efficient visualization via a well-placed set of invariant surfaces. This is because streamlines launched from any smooth curve form a streamsurface by definition. As a consequence, infinitely many streamsurfaces cross through any point of the flow domain. It is thus unclear which (if any) of these surfaces should be selected even locally as a representative of the vector field topology. As a consequence, flow visualization packages tend to rely on user input for seed points of streamlines and streamsurfaces.
Most related ongoing research in the scientific visualization community focuses either on the more accurate computation of streamsurfaces from select seed points (Hultquist Reference Hultquist1992) or streamline selection for streamsurfaces based on visual optimization (Born et al. Reference Born, Wiebel, Friedrich, Scheuermann and Bartz2010; Schulze et al. Reference Schulze, Martinez-Esturo, Günther, Róssl, Seidel, Weinkauf and Theisel2014). A relatively recent realization is that the streamsurfaces framing the instantaneous flow behaviour most efficiently are the key invariant manifolds of the instantaneous velocity field. Local stable and unstable manifolds near stagnation points and closed streamlines indeed well illustrate the instantaneous local velocity field geometry (see, e.g. Peikert & Sadlo Reference Peikert and Sadlo2009) but generally stretch and fold globally. These globally filamented streamsurfaces then lose their ability to demarcate different flow regions efficiently, with unavoidable inaccuracies also arising in their computation (Sadlo & Peikert Reference Sadlo and Peikert2007). Closest in spirit to our work is the observation of Van Wijk (Reference Van Wijk1993), who seeks streamsurfaces as level sets of a scalar function. After defining a scalar distribution within an inflow boundary, the level curves of this scalar field are propagated along streamlines into the flow. As a consequence, the resulting surfaces will generally stretch, fold and self-accumulate, resulting in filamenting streamsurfaces that also depend on the choice of the initial scalar distribution. Reviews of all these approaches in the scientific visualization community can be found, e.g. in Peikert & Sadlo (Reference Peikert and Sadlo2009) and Martinez–Esturo et al. (Reference Martinez–Esturo, Schulze, Rössl and Theisel2013).
A representative set of streamsurfaces for flow illustration should arguably include at least one surface for each typically observed streamline topology, as well as the (generally unobserved) surfaces separating these different topological classes. Finding such a set of streamlines is straightforward for 2-D incompressible flows by the existence of a first integral (conserved quantity) for the equation of instantaneous streamlines. This first integral is the streamfunction, whose set of level curves contains all typically observed streamline families as well as separatrices among them. Therefore, one can systematically scan through the one-parameter family of level curves and select those that stand out to be included in the visualization.
For the lack of a streamfunction in 3-D flows, the above program can only be carried out for streamsurfaces of integrable 3-D flows, such as steady Euler flows that do not satisfy the Beltrami property. For such flows, the Bernoulli function provides a non-degenerate first integral (Arnold & Khesin Reference Arnold and Khesin1999) whose one-parameter family of level surfaces can be scanned and filtered to obtain the required representative set of streamsurfaces. A similar result is available for incompressible velocity fields with a volume-preserving symmetry group (Haller & Mezić Reference Haller and Mezić1998). For general 3-D incompressible flows, however, no first integrals will exist. This is because even steady 3-D flows can be chaotic (see, e.g. Dombre et al. Reference Dombre, Frisch, Greene, Hénon, Mehr and Soward1986) and hence cannot have smooth non-trivial conserved quantities. This is equally valid for the autonomous 3-D differential equation defining the instantaneous streamlines of a generic 3-D unsteady flow, excluding the existence of a regular foliation of the flow domain by streamsurfaces that are level sets of a smooth scalar function.
Often, however, the most important parts of a flow turn out to be vortical regions filled with tubular or toroidal streamsurfaces. One cannot expect these surfaces to necessarily form a continuous family, especially the toroidal ones. Indeed, families of 2-D tori composed of streamlines typically form Cantor sets (as opposed to a continuous family) in the 3-D set of differential equations generating the streamlines (Cheng & Sun Reference Cheng and Sun1989). While this is also the case for classic Hamiltonian systems admitting families of invariant tori (Arnold Reference Arnold1989), those systems nevertheless turn out to be integrable restricted to this Cantor set of tori in an appropriate sense (see Chierchia & Gallavotti Reference Chierchia and Gallavotti1982; Pöschel Reference Pöschel1982). Specifically, smooth functions can be constructed that act as first integrals over the Cantor family of tori but not in the small gaps among those tori.
Motivated by this idea of integrability over Cantor sets in Hamiltonian systems, we seek here smooth scalar functions that serve as approximate first integrals over a set of streamsurfaces forming vortical (elliptic) regions of a given vector field. As there is no widely accepted definition of a vortex, we use the term ‘vortical’ loosely to describe families of toroidal or cylindrical surfaces to which either the velocity, vorticity or barrier field is tangent. The approximate first integrals arising from this procedure will be steady for steady vector fields and time-varying for unsteady vector fields. We construct these (approximate) first integrals by seeking scalar functions whose gradient vector field is as close to being normal to the given vector field as possible. To avoid the trivial solution to this problem, we use a constrained minimization approach that does not allow for globally constant first integrals.
Our method resembles that of Yang & Pullin (Reference Yang and Pullin2010) for the construction of VSFs in inviscid, analytic and highly symmetric flows. Ours, however, differs crucially in that we work with a grid in the physical space over which we expand the unknown approximate first integral in a Fourier series. We then only use the known values of the vector field at the gridpoints. We find that this approach results in an homogeneous linear system of equations whose unique, unit-norm least-squares solution yields the unknown Fourier coefficients of the approximate first integral. Thus, in contrast to Yang & Pullin (Reference Yang and Pullin2010), this procedure is free from any symmetry assumptions on the first integral and does not require rewriting the homogeneous system as an inhomogeneous one under further assumptions. As we show in one of the appendices, these features of our approach significantly enhance the quality of the final solution.
Outside elliptic regions (i.e. in hyperbolic streamline domains), the streamlines are generically chaotic and hence will admit only trivial approximate first integrals. Accordingly, we expect the approximate first integrals to be nearly constant outside vortical regions, while admitting non-trivial shapes inside such regions. In those elliptic regions, level sets of the approximate first integrals will be close to streamsurfaces that form vortical features.
The result of this approach is an automated numerical visualization method that does not require the user to guess seed points in the flow for streamsurfaces in vortical regions. The simplicity and generality of the proposed method allows us to employ it to complex flows defined either analytically or through numerical data. Our examples include spatially periodic integrable and non-integrable flows, a non-periodic vortex ring flow, a V-junction flow and a fully developed turbulent channel flow. The latter flow exemplifies a case wherein exact streamsurfaces tend to obscure the visualization of the most prominent features of the barrier field. Indeed, for such a flow, vortical regions have been delineated via diagnostic tools such as the active version of the finite-time Lyapunov exponent (aFTLE), as described by Haller et al. (Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020). Yet, tracking barrier streamlines originating in the neighbourhood of aFTLE ridges quickly results in the streamsurface falling apart, despite some initial vortical motion. In contrast, the structures based on the approximate first integral are able to follow closely the valleys around the aFTLE ridges allowing for a better visualization.
To obtain these results, we divide the computational domain into smaller subdomains and seek approximate first integrals in each one of them separately. In this way, we can capture the most salient structures up to a given scale without the problem of obstruction by smaller structures. These in turn can be captured by further refining the domain subdivision based on their signatures in the aFTLE field. Further, upon assuming mild convexity for the vortical structures to be extracted, we obtain families of barrier surfaces whose outermost members act as vortex boundaries. Thus, we propound a way out of the isocontour value dilemma that besets the typical vortex identification criteria used in the literature (Kim, Moin & Moser Reference Kim, Moin and Moser1987; Jeong & Hussain Reference Jeong and Hussain1995; Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999), causing them to produce dissimilar structures for different isocontour values.
2. Set-up of the minimization scheme and reconstruction algorithm
Let 
$\boldsymbol {v}(\boldsymbol {x},t)$ be a smooth vector field defined on some open subset 
$U\subseteq \mathbb {R}^{3}$. The associated dynamical system for the instantaneous streamlines of this vector field at time 
$t$ is given by
\begin{equation} {\boldsymbol{x}}^\prime = \boldsymbol{v}(\boldsymbol{x},t),\quad \boldsymbol{x} \in U, \end{equation}
where the prime denotes differentiation with respect to a curve-parameter 
$s\in {\mathbb {R}}$ along the streamline. A continuously differentiable, scalar function 
$H(\boldsymbol {x},t)$ is called an (instantaneous) first integral for 
$\boldsymbol {v}(\boldsymbol {x},t)$ at time 
$t$ if it is constant along each solution of (2.1), i.e. 
$({\partial }/{\partial s}) H(\boldsymbol {x}(s),t) = 0$. This condition implies that
\begin{equation} \boldsymbol{\nabla} H(\boldsymbol{x},t)\boldsymbol{\cdot} \boldsymbol{v}(\boldsymbol{x},t)=0, \end{equation}
for all 
$\boldsymbol {x}\in U$, where 
$\boldsymbol {\nabla }$ denotes the gradient with respect to 
$\boldsymbol {x}$.
For an arbitrary vector field 
$\boldsymbol {v}$, no first integral will exist, in general, pointwise, i.e. for all 
$\boldsymbol {x}\in U$. We may, however, relax the constraint (2.2) by seeking a function 
$H(\boldsymbol {x},t)$ that minimizes the functional
\begin{equation} J[H]=\frac{1}{2}\int_{U}|\boldsymbol{\nabla} H \boldsymbol{\cdot} \boldsymbol{v}|^{2}\,{\rm d}V, \end{equation}
which measures the average deviation of 
$H(\boldsymbol {x},t)$ from being an exact, pointwise first integral in the domain 
$U$ at time 
$t$. Any minimizer of (2.3) is called an approximate first integral.
First, let us assume that the domain 
$U$ is triply periodic. This allows us to expand 
$H$ in a Fourier series and define its modal truncation of order 
$N$ as
\begin{equation} H^{{\leq} N}(\boldsymbol{x})=\sum_{0<|\boldsymbol{k}|\leq N}\hat{H}_{\boldsymbol{k}}{\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}. \end{equation}For the application to the examples in later sections, let us remark that the truncated expansion (2.4) can also be used locally on non-periodic domains, as long as we stay away from the domain boundaries. The integrand of the functional (2.3) in Fourier basis takes the form
\begin{equation} \boldsymbol{\nabla} H\boldsymbol{\cdot}\boldsymbol{v}=\sum_{0<|\boldsymbol{k}|\leq N}\hat{H}_{\boldsymbol{k}} {\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}. \end{equation}Because of the linearity of the gradient operator, the expression in (2.5) is linear in the unknown Fourier coefficients.
We assume that the vector field 
$\boldsymbol {v}$ is known on a discrete, three-dimensional grid of points which we enumerate from 
$1$ through 
$m$, where 
$m$ is the total number of points. Also, to target data-driven applications specifically, we work with the discretized version of (2.3) using the vector of pointwise inner products defined as
\begin{equation} \left[\begin{array}{cccc} \boldsymbol{\nabla} H_{1}\boldsymbol{\cdot}\boldsymbol{v}_{1} & \boldsymbol{\nabla} H_{2}\boldsymbol{\cdot}\boldsymbol{v}_{2} & \cdots & \boldsymbol{\nabla} H_{m}\boldsymbol{\cdot}\boldsymbol{v}_{m}\end{array}\right]^{{\rm T}}=\boldsymbol{\mathsf{C}} \boldsymbol{\boldsymbol{h}}, \end{equation}
with 
$\boldsymbol{\mathsf{C}}\in {\mathbb {C}}^{m\times n}$. Here, 
$n$ is the number of modes used, 
${\mathsf{C}}_{ij}=\exp ({{\rm i}\boldsymbol {k}_{j}\boldsymbol {\cdot}\boldsymbol {x}_{i}})\boldsymbol {k}_{j}\boldsymbol {\cdot}\boldsymbol {v}_{i}$ is the 
$(i,j)$ entry of 
$\boldsymbol{\mathsf{C}}$ and 
$\boldsymbol {\boldsymbol {h}}=\{ \hat {H}_{\boldsymbol {k}}|\boldsymbol {k}\in \mathbb {Z}^{3},\, 0<|\boldsymbol {k}|\leq N\}$ is the vector comprising the Fourier coefficients to be determined. Approximating a first integral through the functional (2.3) then amounts to minimizing the squared vector norm 
$|\boldsymbol{\mathsf{C}}\boldsymbol {\boldsymbol {h}}|^2$. To exclude the trivial solution 
$H\equiv {\rm const.}$ from our analysis, we add the constraint 
$|\boldsymbol {\boldsymbol {h}}|^2=1$.
Solving this optimization problem is equivalent to finding the eigenvector corresponding to the minimal eigenvalue of the symmetric matrix 
$\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{C}}^{*}\boldsymbol{\mathsf{C}}$, where 
$\boldsymbol{\mathsf{C}}^*$ denotes the conjugate transpose of 
$\boldsymbol{\mathsf{C}}$ (see Appendix B). Since the eigenvectors of 
$\boldsymbol{\mathsf{A}}$ are the right-singular vectors of 
$\boldsymbol{\mathsf{C}}$, the solution to our algorithm can also be calculated from the singular-value decomposition (SVD) of 
$\boldsymbol{\mathsf{C}}$.
We refrain from expanding the known vector field 
$\boldsymbol {v}$ into a Fourier series and we write down (2.2) explicitly for every point in the computational grid without working with the coefficients of each Fourier mode. These features distinguish our method from the one presented by Yang & Pullin (Reference Yang and Pullin2010), as already noted in § 1. This difference will allow us to obtain unique solutions as well as apply our approach even to turbulent flows without additional assumptions.
In the absence of further constraints, the resulting approximate first integral 
$H$ will, generally, be a complex-valued scalar field. Denoting by 
$H_{r}$ and 
$H_{i}$ its real and imaginary parts, respectively, we have 
$\langle \boldsymbol {\nabla } H,\boldsymbol {v} \rangle _{\mathbb {C}} = \langle \boldsymbol {\nabla } H_{{r}},\boldsymbol {v} \rangle - \langle \boldsymbol {\nabla } H_{{\rm i}},\boldsymbol {v} \rangle \mkern 3mu {\rm i}$. This implies that by considering the expression in (2.6), we essentially optimize both the real and imaginary parts of 
$H$. This allows us to use 
$| H |$ as the approximate first integral. Alternatively, we can require 
$H$ to be real a priori. We discuss the latter procedure in Appendix C and show that its results are similar to those obtained without imposing this constraint.
As already noted, for generic flows, a non-trivial approximate first integral is only expected to exist in vortical regions. Outside such regions, we expect our algorithm to yield almost flat first integrals. Such a function, locally supported on several vortical regions, would generally require a large number of Fourier basis functions, which in turn would lead to numerical inefficiencies and cost. To avoid this, we work with a comparably small number of basis functions and only use the level surfaces of the emerging approximate first integral in regions where those surfaces are indeed nearly tangent to the given vector field. To identify such regions, we introduce the invariance error as
\begin{equation} E_{A} =\frac{1}{A} \sum_{i=1}^{p} \left| \frac{\boldsymbol{\nabla} \left| H_{i} \right| \boldsymbol{\cdot} \boldsymbol{v}_{i}}{\left|\boldsymbol{\nabla} \left|H_{i}\right| \right| \left| \boldsymbol{v}_{i} \right|} \right|, \end{equation}
where 
$A$ is the surface area of a level set and 
$p$ the number of points in the level set. This type of error estimate was first introduced by Yang & Pullin (Reference Yang and Pullin2010). In our visualizations, we will discard streamsurface candidates with invariance errors exceeding a certain threshold value.
Finally, we note that our minimization procedure can also be viewed as finding for 
$\boldsymbol {v}$, in the appropriate norm, the closest member of the integrable, 3-D incompressible vector field family
\begin{equation} {\boldsymbol{x}}^\prime = \boldsymbol{\mathsf{J}}(\boldsymbol{x})\boldsymbol{\nabla} H(\boldsymbol{x},t),\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{J}}\equiv\boldsymbol{0}, \end{equation}
with 
$\boldsymbol{\mathsf{J}}=-\boldsymbol{\mathsf{J}}^{\rm T}$ and 
$\boldsymbol {\nabla }\boldsymbol {\cdot} \boldsymbol{\mathsf{J}}$ denoting the divergence of the tensor field 
$\boldsymbol{\mathsf{J}}$. Indeed, all these vector fields in (2.8) share the same streamsurfaces, the level sets of 
$H$. Working with (2.8) directly, however, is much more demanding numerically in our experience and would also require specific assumptions on the form of 
$\boldsymbol{\mathsf{J}}$.
Before moving to specific examples, we note that finding an exact, pointwise first integral in (2.2) is not a well-posed problem by itself. Indeed, if 
$H(\boldsymbol {x},t)$ is a solution, then, for any sufficiently smooth function 
$F$, 
$F(H(\boldsymbol {x},t))$ will also be a solution due to the homogeneity of (2.2). This would not be an issue for the detection of streamsurfaces if the isocontours of 
$H$ and 
$F(H)$, expressed through the gradient of these fields, represent the same geometric and topological features for the given vector field 
$\boldsymbol {v}$. Unfortunately, we can construct simple counter-examples where this is not the case. For instance, if we denote by 
$v_x, v_y$ and 
$v_z$ the three components of 
$\boldsymbol {v}$ and assume that 
$v_x = 0$, then the function 
$G(x) H(y,z)$ will be an exact, pointwise first integral as long as 
$v_y \theta _y H + v_z \theta _z H = 0$. By tweaking 
$G$, one can easily obtain markedly different topological features resulting from the corresponding streamlines. We also refer to Pullin & Yang (Reference Pullin and Yang2014) for more examples of first integrals with different topology for the vorticity field of the Taylor–Green flow stemming from the superposition of independent solutions to (2.2).
Consequently, our variational principle (2.3) will exhibit the same non-uniqueness issues whenever an exact, pointwise first integral is admitted by the underlying vector field. To resolve this, we will only consider approximate first integrals for which the weakest eigenvalue of 
$\boldsymbol{\mathsf{A}}=\boldsymbol{\mathsf{C}}^*\boldsymbol{\mathsf{C}}$ is not considered (numerically) zero. In addition, we will only retain the streamsurfaces whose topology remains unaltered when, for the same set of grid points, a larger number of Fourier modes is used allowing for small geometric refinements due to the increased accuracy of the solution. If these two conditions are met, we will consider the resulting structures as robust and they will be included in the visualization. Furthermore, we will see that the more complex a flow is, the larger the spectral gap to the second-smallest eigenvalue will be. Irrespective of this gap, however, in all the examples that follow, we will build our solution based only on the eigenvector corresponding to the smallest eigenvalue of 
$\boldsymbol{\mathsf{A}}$. We close this section by noting that our approach is in agreement with the findings of Xiong & Yang (Reference Xiong and Yang2017, Reference Xiong and Yang2019) who showed that the construction of VSFs in turbulent flows through the use of partial differential equations leads to robust structures despite the non-uniqueness issues.
3. Approximate first integrals for explicit solutions of the Euler equations
In this section, we illustrate our minimization algorithm to construct approximate first integrals and use their level sets as approximate streamsurfaces in analytic examples.
3.1. ABC flow
As a first test case, we investigate the ABC (Arnold–Beltrami–Childress) class of flows (Dombre et al. Reference Dombre, Frisch, Greene, Hénon, Mehr and Soward1986; Henon Reference Henon1966), defined as
\begin{equation} \left. \begin{gathered} \displaystyle \dot{x} = A \sin z + C \cos y, \\ \displaystyle \dot{y} = B \sin x + A \cos z, \\ \displaystyle \dot{z} = C \sin y + B \cos x, \end{gathered} \right\} \end{equation}
for 
$A,B,C\in \mathbb {R}$ and 
$(x,y,z) \in [0,2{\rm \pi} ]^3$, together with periodic boundary conditions. The right-hand side of (3.1) defines an exact steady solution to the incompressible Euler equations. For 
$ABC \neq 0$, the flow exhibits chaotic behaviour (Dombre et al. Reference Dombre, Frisch, Greene, Hénon, Mehr and Soward1986; Henon Reference Henon1966), whereas some analytic non-integrability results have been reported by Ziglin (Reference Ziglin1988, Reference Ziglin1998).
3.1.1. Integrable case
We first analyse the ABC flow with 
$A=0$ for which (3.1) is completely integrable. For 
$BC \neq 0$, an exact, pointwise first integral is given by 
$H_1(x,y) = C \sin y + B \cos x$, while another independent first integral can be constructed through the use of elliptic functions (Llibre & Valls Reference Llibre and Valls2012). The level sets of 
$H_1$ are depicted in figure 1(a) on one cross-section of the computational domain for 
$B=\sqrt {2}$ and 
$C = 1$. These curves, therefore, represent the intersections of a representative set of streamsurfaces with the 
$z=0$ plane. These streamsurfaces are also exact invariant manifolds for the Lagrangian particle motions in this steady flow.
Analysis of the integrable ABC flow using a computational grid of 
$100^3$ points and approximately 
$9000$ Fourier modes. (a) Intersections of the level surfaces of 
$H_1$ with the 
$z=0$ plane. (b) Intersections of the level surfaces of the approximate first integral 
$H$ with the 
$z=0$ plane. (c) Same as panel (b) but after the removal of small-scale structures of panel (b) as well as the structures with 
$E_{l} > 10^{-5}$.
To test our proposed algorithm for constructing approximate first integrals, we start from a discretized version of the full 3-D velocity field (3.1) along 
$100$ points per direction, using approximately 
$9000$ Fourier modes in (2.4)–(2.6). Solving the underlying optimization problem using the algorithm in Appendix B, we obtain the results shown in figure 1(b) at the same cross-section as in figure 1(a). The numerically constructed approximate level sets match the analytic first integral closely. To measure the proximity of streamsurfaces and approximate streamsurfaces along the 
$z=0$ plane, we introduce a planar version of the general invariance error (2.7) by defining
\begin{equation} E_{l} =\frac{1}{l} \sum_{j=1}^{p} \left| \frac{\boldsymbol{\nabla} \left| H_{j} \right| \boldsymbol{\cdot} \boldsymbol{\nabla} H_{1j}}{\left| \boldsymbol{\nabla} \left|H_{j}\right| \right| \left| \boldsymbol{\nabla} H_{1j} \right|} \right|, \end{equation}
where 
$l$ is the length of the level and 
$p$ is the number of points for each level set. We use this metric to remove level curves with fewer than 
$30$ points and those with invariance errors 
$E_{l} > 10^{-5}$. The results shown in figure 1(c) confirm that choosing these thresholds removes small-scale artefacts arising from numerical inaccuracies.
We perform the same analysis on the 
$z=0$ plane using 
$150$ points per direction but only approximately 
$1200$ Fourier modes. The results are depicted in figure 2. Again, we observe close agreement between the known first integral and the reconstructed curves. At some points, the agreement is even closer when compared to figure 1(c) due to the higher spatial resolution, even though the number of Fourier modes is significantly smaller.
Same as figure 1 but with a computational grid of 
$150^3$ points and 
$1200$ Fourier modes.
3.1.2. Non-integrable case
For a different set of parameter values (
$A = \sqrt {3}$, 
$B = \sqrt {2}$ and 
$C = 1$), the ABC flow (3.1) is non-integrable and shows chaotic behaviour in some regions. The dynamic behaviour of trajectories for this set of parameter values is well studied, including the KAM-type tori highlighted by Poincaré maps (Dombre et al. Reference Dombre, Frisch, Greene, Hénon, Mehr and Soward1986) and elliptic Lagrangian coherent structure (LCS) techniques (Blazevski & Haller Reference Blazevski and Haller2014; Oettinger, Blazevski & Haller Reference Oettinger, Blazevski and Haller2016). We use this velocity field as a benchmark to test different solution algorithms for finding an approximate first integral for a non-integrable flow. Also, this will serve as a proof of concept for finding elliptical regions.
We have already noted that the unit-norm least-squares solution to the homogeneous system of (2.6) coincides with the right-singular vector of 
$\boldsymbol{\mathsf{C}}$ associated with its smallest singular value or, equivalently, with the eigenvector associated with the smallest eigenvalue of 
$\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{C}}^{*}\boldsymbol{\mathsf{C}}$. To improve numerical stability, the SVD-based solution is preferred (Golub & Pereyra Reference Golub and Pereyra1973). Indeed, the eigenvalue calculation requires a matrix multiplication to form 
$\boldsymbol{\mathsf{A}}$, which invariably squares the condition number of 
$\boldsymbol{\mathsf{C}}$. For comparison, we calculate both the singular-vector-based and the eigenvector-based solutions on the triply periodic box 
$[0,2 {\rm \pi}]^3$ with 
$100$ points per direction and 
$N=13$ (or 
$9170$ Fourier modes). In this setting, running the SVD algorithm of MATLAB on 
$\boldsymbol{\mathsf{C}}$, which is a 
$(10^2)^3 \times 9170$ matrix, would require an exorbitant amount of memory (more than 
$400$ GB), indicating that the classical SVD algorithm is not optimized for tall-skinny matrices (such as our coefficient matrix). To proceed with our comparison test, we instead follow the modified SVD method discussed in Appendix D for tall-skinny matrices.
With this modification to the SVD-based solution, the results from the two approaches for the non-integrable ABC flow are presented in figure 3 on the 
$z=0$, 
$y=0$ and 
$x=0$ planes. We observe that the differences between the eigenvector-based and singular-vector-based computations are marginal, indicating that the larger condition number of 
$\boldsymbol{\mathsf{A}}$ does not affect the results. Furthermore, we use the same three planes as Poincaré sections to integrate trajectories up to an arclength of 
$10^4$ from a uniform grid of 
$20 \times 20$ initial conditions on each section. We observe a very good agreement between the predicted structures and the intersections of the KAM surfaces with each of these sections. This is highlighted perhaps even better by the reconstructed KAM surfaces as approximate streamsurfaces in figure 4, which are to be contrasted with the 
$\lambda _2$-based structures in Appendix A. Since the ABC flow is a Beltrami flow, its velocity 
$\boldsymbol {v}$ is parallel to the vorticity 
$\boldsymbol {\omega } = \boldsymbol {\boldsymbol {\nabla }} \times \boldsymbol {v}$ (in fact, 
$\boldsymbol {v}=\boldsymbol {\omega }$); consequently, the approximate-first-integral-based tori we have constructed are also VSFs. This illustrates that our algorithm can identify VSFs in flows where the methodology of Yang & Pullin (Reference Yang and Pullin2010) is inapplicable. Indeed, as already noted, the non-integrable ABC flow has chaotic streamlines and, thus, no symmetry assumptions regarding these streamlines can be used to accelerate the convergence rate for the optimization technique presented by Yang & Pullin (Reference Yang and Pullin2010). Even if this rate was irrelevant, however, expanding the known velocity field in a Fourier series would result in an optimization problem with many (numerically) zero eigenvalues and, thus, infinitely many possible minimizers.
Analysis of the non-integrable ABC flow using a computational grid of 
$100^3$ points and 
$9170$ Fourier modes. Level sets of the approximate first integral at (a,d) 
$z=0$, (b,e) 
$y=0$ and (c,f) 
$x=0$. Panels (a–c) are constructed from the eigenvector of 
$\boldsymbol{\mathsf{A}}$ corresponding to the smallest eigenvalue, whereas panels (d–f) are produced using the SVD of 
$\boldsymbol{\mathsf{C}}$. The overlaid Poincaré map (black dots) on each section is based on a uniform grid of 
$20 \times 20$ initial conditions.
Two different views of the approximate streamsurfaces (level sets of the approximate first integral) closely approximate the KAM-type surfaces of the non-integrable ABC flow in elliptic regions. Also shown are iterations of the Poincaré map (black dots) on three orthogonal planes. The results were obtained using the weakest eigenvector of the positive definite matrix 
$\boldsymbol{\mathsf{A}}$. See also the supplementary movie 1 available at https://doi.org/10.1017/jfm.2022.992.
Upon taking a closer look at the results of figure 3, we notice that the reconstructed level sets attain their values in a longer range (i.e. 
$[0,3.5]$) for the larger KAM surfaces, whereas, in the vicinity of the smaller structures, they are confined to a narrow band (i.e. 
$[0.25,0.35]$). Here the adjectives larger and smaller are used to refer to either the area (figure 3) or the volume (figure 4) these families enclose. This is a type of overfitting that we would like to mitigate. One way to achieve this is by considering a slightly different optimization problem (see Appendix E) which resembles the one put forward by Yang & Pullin (Reference Yang and Pullin2010). This different approach, however, turns out to be computationally intense, posing severe limitations on its use for typical grid sizes while its results are arguably of inferior quality. All in all, obtaining the solution to the proposed algorithm as the eigenvector corresponding to the smallest eigenvalue of 
$\boldsymbol{\mathsf{A}}$ is computationally superior to all the other techniques used and, thus, it is the one that we will follow in the rest of this article.
We conclude this section by performing a convergence analysis for different numbers of modes in figure 5. We observe that the least-squares error (as the smallest eigenvalue of 
$\boldsymbol{\mathsf{A}}$) approaches zero as the number of modes used in the analysis increases. Furthermore, the second smallest eigenvalue converges to the smallest one for higher modes. In Appendix B, we show that when the smallest eigenvalues are almost equal, we can construct the solution as a linear combination of the eigenvectors corresponding to these near-identical eigenvalues. Here, however, we can base our solution only on the weakest eigenvector since we have 
$d_1 = 0.015415$ and 
$d_2=0.01596$ for 
$N=13$. This gap between 
$d_1$ and 
$d_2$ will prove significantly larger in the following, data-based examples, confirming the uniqueness of solutions to the optimization problem.
Five smallest eigenvalues of 
$\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{C}}^{*}\boldsymbol{\mathsf{C}}$ for different numbers of modes (
$2108, 3070, 4168, 5574, 7152, 9170$ modes for 
$N=8,9,10,11,12,13$, respectively).
3.2. Further analytic solution to the Euler equations
A set of analytic, unsteady, tri-periodic laminar solutions of the Navier–Stokes equations was put forward by Antuono (Reference Antuono2020). Here, we consider only the steady part of these solutions, which is a Beltrami solution to the Euler equations with no known first integral. We therefore expect the streamlines of this velocity field to be chaotic and the overall dynamics to be non-integrable. Representative (approximate) streamsurfaces have not yet been constructed for this flow in the literature.
The velocity field is given by
\begin{equation} \boldsymbol{v}=\dfrac{4\sqrt{2}}{3\sqrt{3}}\left(\begin{array}{c} \sin\left(x-\dfrac{5{\rm \pi}}{6}\right)\cos\left(y-\dfrac{\rm \pi}{6}\right)\sin(z)-\cos\left(z-\dfrac{5{\rm \pi}}{6}\right)\sin\left(x-\dfrac{\rm \pi}{6}\right)\sin(y)\\ \sin\left(y-\dfrac{5{\rm \pi}}{6}\right)\cos\left(z-\dfrac{\rm \pi}{6}\right)\sin(x)-\cos\left(x-\dfrac{5{\rm \pi}}{6}\right)\sin\left(y-\dfrac{\rm \pi}{6}\right)\sin(z)\\ \sin\left(z-\dfrac{5{\rm \pi}}{6}\right)\cos\left(x-\dfrac{\rm \pi}{6}\right)\sin(y)-\cos\left(y-\dfrac{5{\rm \pi}}{6}\right)\sin\left(z-\dfrac{\rm \pi}{6}\right)\sin(x) \end{array}\right). \end{equation}
Using a uniform grid of 
$20 \times 20$ initial conditions on the 
$y=0$ plane, we integrate trajectories up to an arclength of 
$10^4$ and, upon retaining their long-term behaviour from the interval 
$[5 \times 10^3, 10^4]$, we show the resulting Poincaré map in figure 6(a) where a plethora of KAM tori are discernible. Selecting the triply periodic box 
$[0,2 {\rm \pi}]^3$ and sampling it with 
$100^3$ points, we run our algorithm for 
$N=13$. Upon constructing isosurfaces for 
$10$ different isovalues of the resulting approximate first integral, we locate their intersections with the 
$y=0$ plane and superimpose them on figure 6(a). We then compute the invariance error based on (2.7) and colour the isocontours of figure 6(a) blue or red depending on whether their average error corresponds to an angle of less or more than 
$5^{\circ }$. Based on this, contours lying inside the chaotic sea of the Poincaré map show the largest invariance error despite our algorithm not using any knowledge of the long-term, chaotic dynamics.
(a) Comparison of the Poincaré map on the plane 
$y=0$ for the steady Euler flow (3.3) overlaid on the intersections of the tori obtained from an approximate first integral with the same plane for 
$N = 13$. The blue (red) isocontours depict tori whose invariance error (see (2.7)) is smaller (larger) than 
$5^{\circ }$ on average. (b) 95th percentile of the distance (in non-dimensional units) between trajectories emanating from the isocontours of panel (a) with the corresponding 2-D tori inside the computational box 
$[0,2{\rm \pi} ]^3$. The points to the left (right) of the dash–dotted line correspond to trajectories originating in the blue (red) contours of panel (a).
For every isocontour of figure 6(a), we launch trajectories following the vector field 
$\boldsymbol {v}$ of (3.3) until they leave the computational domain 
$[0,2{\rm \pi} ]^3$, and then calculate their distances from the corresponding isosurfaces. Thus, in figure 6(b), we present the 95th percentile of those distances for each streamsurface after separating them according to figure 6(a), i.e. we depict the trajectories corresponding to the blue contours in the left side of the dash–dotted line, whereas the ones for the red contours in the right side. We note the correlation between the retained (blue) isocontours of figure 6(a) and the significantly smaller percentiles of figure 6(b).
Similarly, we present the reconstructed tori for 
$N = 15$ or 
$14\,146$ modes in figure 7(a). We observe that for 
$N=13$, the approximate first integral captures virtually all the KAM surfaces indicated by the Poincaré map. In contrast, for 
$N=15$, some of the structures are captured more accurately while others are missed completely. The convergence analysis depicted in figure 8(a) shows that the least-squares error follows a declining trend as 
$N$ grows. This prompts us to consider an error measure similar to the one in (2.7), defined as
\begin{equation} E_m=\frac{1}{m}\sum_{i=1}^{m}\left| \frac{\boldsymbol{\nabla} \left| H_{i} \right| \boldsymbol{\cdot} \boldsymbol{v}_{i}}{\left| \boldsymbol{\nabla} \left| H_{i} \right|\right| \left| \boldsymbol{v}_{i}\right| } \right|. \end{equation}
In this expression, the summation is taken over all the grid points, providing an estimate for the mean invariance error of the entire solution. This allows us to make a direct comparison among solutions corresponding to different numbers of modes. Specifically, excluding points that lie in the vicinity of fixed points of either 
$\boldsymbol {v}$ or 
$\boldsymbol {\nabla } | H |$, we show the dependence of the invariance error 
$E_m$ on the number of Fourier modes used in our algorithm in figure 8(b). We observe that the error attains the minimum value for 
$N=13$, in agreement with what is inferred from figure 6(a).
(a) Same as in figure 6(a) with the tori reconstructed for 
$N = 15$. (b) Closeup view on a region filled with two families of invariant tori. Overlaid on the Poincaré map are the tori obtained from an approximate first integral for 
$N = 19$.
Moreover, to mitigate minor discrepancies between the reconstructed tori and those outlined by the Poincaré map of figure 6(a), we increase the number of Fourier modes to 
$28\,670$ (
$N=19$), while keeping the same set of grid points. A closeup view of a region filled with invariant tori in figure 7(b) confirms that there is a close agreement with the tori obtained from an approximate first integral. There is, however, a trade-off between the increased accuracy and the ensuing computational burden that an end user must consider.
We also note that, while for 
$N=15$ we have 
$d_1 \approx d_2 = 839.8079$, for 
$N=19$ we have 
$d_1 = 204.5165$ and 
$d_2 = 218.2078$, further corroborating that our approach leads to a unique solution. Finally, we conclude this section with figure 9(a) showing 3-D rendered images of the three Poincaré maps on the planes 
$x=y=z=0$, figure 9(b) showing representative streamsurfaces obtained as level surfaces of an approximate first integral for 
$N=13$ and figure 9(c) showing the superimposition of panel (a) on panel (b), which confirms the close agreement between the expected and reconstructed structures.
3.3. Hill's spherical vortex
We now turn to a spatially non-periodic, integrable flow given by Hill's spherical vortex (Hill Reference Hill1894). The axial-symmetric (approximately the 
$z$-axis) streamfunction for Hill's solution to the Euler equations is
\begin{equation} \psi(h,z)=\left\{ \begin{array}{@{}ll} \dfrac{3}{4}U_{0}h^{2}\left(1-r^{2}\right), & r\leq1,\\ -\dfrac{1}{2}U_{0}h^{2}\left(1-\dfrac{1}{r^{3}}\right), & r>1, \end{array} \right. \end{equation}
where 
$h$ is the distance from the axis of symmetry (
$h^2 = x^2 +y^2$) and 
$r^{2}=h^{2}+z^{2}$. Using the relations 
$u_{h}=-({1}/{h})({\partial \psi }/{\partial z})$ and 
$u_{z}=({1}/{h})({\partial \psi }/{\partial h})$, we obtain the corresponding velocity field in Cartesian coordinates
\begin{equation} \displaystyle \boldsymbol{v}(x,y,z)=\left\{ \begin{array}{@{}ll} \dfrac{3}{2}U_{0} \left(xz, \ yz, \ 1-(r^2+h^2)\right) , & r\leq1,\\ \dfrac{3}{2}U_{0} \left(\dfrac{xz}{r^5}, \ \dfrac{yz}{r^5}, \ \dfrac{2}{3} \dfrac{1}{r^3} - \dfrac{h^2}{r^5}\right), & r>1. \end{array} \right. \end{equation} As for periodic domains, we can compute a Fourier expansion on a bounded subdomain, bearing in mind that the convergence of the partial Fourier sum will be slow in general (Gottlieb & Shu Reference Gottlieb and Shu1997). To mitigate this issue, we will have to consider a sufficiently high number of modes. Furthermore, due to the well-known Gibbs phenomenon, there will be sizeable spurious oscillations in the approximate first integral near the box boundary that do not diminish after an increase in the number of modes. To damp this effect in our algorithm, we discard all the reconstructed surfaces that fall within 
$5\,\%$ of the box size in all three directions.
The velocity field of (3.6) has toroidal streamsurfaces inside the spherical domain 
$\{\boldsymbol {x} \in \mathbb {R}: |\boldsymbol {x}| \leq 1\}$ (see figure 10a). We reconstruct these streamsurfaces over the computational domain 
$[-2,2]^3$ with 
$60$ points per direction, running our algorithm for different numbers of Fourier modes with 
$U_{0}=1$. In some cases with 
$N<15$, the non-periodic nature of the velocity field results in the breakdown of the reconstructed toroidal structures. For 
$N\ge 15$, however, we obtain fully symmetric solutions. Specifically, for illustration purposes, we work with 
$N=17$ (or 
$20\,478$ Fourier modes) which yields a unique solution corresponding to the smallest eigenvalue of 
$\boldsymbol{\mathsf{A}}$, 
$d_{1} = 23.5169$ (
$d_2 = 25.3118$). Based on this, we create 
$10$ isosurfaces in the interval 
$\big[ |H|_{min},|H|_{max} \big]$. Upon locating the intersections of these surfaces with eight radially equidistant planes, we launch trajectories corresponding to these intersections and compute the pointwise distance of the solution curves to the reconstructed surfaces. The results (in terms of percentiles) are presented in figure 10(b). A comparison between the reconstructed streamsurfaces and indicative solution curves is given in figure 11.
(a) Level sets of Hill's streamfunction on the 
$x=0$ plane. (b) 25th, 50th, 75th and 95th percentile of the pointwise distances (in non-dimensional units) between solution curves and 
$10$ different reconstructed streamsurfaces for 
$N=17$.
3.4. Flow inside a V-junction
Next we investigate the flow inside a V-junction, as depicted in figure 12. Despite the simple geometry, recent experiments have suggested that pumping a particle-laden fluid into such a configuration allows light particles, such as gas bubbles in water, to be permanently trapped in the junction (Vigolo, Radl & Stone Reference Vigolo, Radl and Stone2014). This phenomenon arises for a wide range of Reynolds numbers and for various junction angles (Ault et al. Reference Ault, Fani, Chen, Shin, Gallaire and Stone2016). Here we consider one such flow with junction angle 
$70^{\circ }$ and 
${Re} = ( \bar {u}L/\nu ) = 230$ (Shin, Ault & Stone Reference Shin, Ault and Stone2015), where 
$\bar {u}$ is the average inlet flow speed, 
$L$ is the side length of the square channel and 
$\nu$ is the kinematic viscosity.
Velocity streamlines, colour coded with their magnitude, emanating from the inlet of the V-junction and leaving the domain from the two outlets. For 
${Re} = 230$ and a junction angle of 
$70^{\circ }$, four symmetric vortex-breakdown bubbles are portrayed in grey.
We use a finite-volume solver from the OpenFOAM library to obtain a steady solution to the 3-D incompressible Navier–Stokes equation (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998). The same numerical solution has recently been analysed using methods from dynamical systems theory, which revealed large, anchor-shaped trapping regions for light particles (Oettinger et al. Reference Oettinger, Ault, Stone and Haller2018). These trapping regions, however, have been invariably linked to bubble-type vortex breakdown structures in the fluid flow, which are formed downstream in the junction (Ault et al. Reference Ault, Fani, Chen, Shin, Gallaire and Stone2016). These structures are depicted as dark grey blobs in figure 12, obtained from a careful advection of streamlines from the vicinity of known stagnation points. Their construction, therefore, is by no means automated and assumes a detailed knowledge of the streamline geometry.
Each vortex breakdown bubble is demarcated by the 2-D stable and unstable manifolds of two saddle-type fixed points. In a generic 3-D flow, these two manifolds do not coincide because that configuration would not be structurally stable. Instead, they are expected to intersect transversely. In such a scenario, the original streamsurface breaks down allowing fluid particles from the main stream to be entrained into the bubble and, conversely, particles to return from the bubble to the main stream (Holmes Reference Holmes1984; Peikert & Sadlo Reference Peikert and Sadlo2007). Sotiropoulos, Ventikos & Lackey (Reference Sotiropoulos, Ventikos and Lackey2001) showed that a very careful mesh refinement is required to reveal the splitting of the manifolds. Here, we will only be interested in reconstructing a closed streamsurface that manifests minimal fluid exchange (see Peikert & Sadlo (Reference Peikert and Sadlo2007) for a discussion of this surface).
We select as computational domain the box depicted in figure 13 with 
$110, 90, 80$ points in the 
$x,y$ and 
$z$ direction, respectively. This box is chosen so that the bubble-like vortical structure lies approximately in its centre. The dimensions of the box are 
$L_{x} = 1\, \textrm {m}$, 
$L_{y} = 0.6\, \textrm {m}$ and 
$L_{z} = 0.4\, \textrm {m}$. This procedure incorporates a priori information approximating the rough location of the streamsurface of interest. We will see in the next section how the same algorithm can be extended to uncover a priori unknown structures in a turbulent flow.
Different views of the computational box used to construct an approximate first integral for one of the bubble-like structures in the V-junction flow.
We use our algorithm with 
$N=15$ or 
$14\,146$ Fourier modes and show a 2-D cross-section at the middle of the computational domain in figure 14(a). For the smallest eigenvalue of 
$\boldsymbol{\mathsf{A}}$, 
$d_1 = 58.1538$ (
$d_2 = 376.9696$), a very pronounced circular structure is clearly visible approximately in the middle of the domain bordering an otherwise flat landscape. To obtain the full, 3-D reconstruction, we use 
$40$ equidistant values between 
$| H |_{min}$ and 
$| H |_{max}$ and extract 
$28$ different level surfaces of the approximate first integral. We then sort them in decreasing order based on the volume they enclose.
Results for an approximate first integral in the V-junction flow, framing one of the elliptical vortex regions. (a) Approximate first integral distribution on a plane which coincides with the middle, in the 
$x$ direction, of the computational domain presented in figure 13. (b) Normalized invariance error as a function of the extracted isosurfaces sorted in descending order with respect to their volume.
Plotting this error estimate for the extracted isosurfaces (figure 14b) yields a global minimum for the surface 
$i=25$ and a local minimum for the surface 
$i=15$. Rendered depictions of the extracted isosurfaces corresponding to these two minima are shown in figure 15. We note the close agreement between the reconstructed structures and the structure delineated by the judiciously chosen streamlines. This is in stark contrast to the structures suggested by the 
$Q$-criterion as demonstrated in Appendix A.
Streamsurfaces as level sets of an approximate first integral in the V-junction flow, corresponding to one local 
$(i = 15)$ and the global minimum 
$(i = 25)$ of the normalized error 
$E_{A}$ depicted in figure 14(b). The duct is cut transversely to help the visualization.
4. Momentum transport barriers in a turbulent channel flow
As we have mentioned in § 1, the shape of streamsurfaces is observer-dependent: their geometry changes under general time-dependent rotations and translations of the observer. Therefore, unless the flow has a distinguished frame in which streamsurfaces coincide with material surfaces, streamsurfaces simply highlight features of the velocity field in a given frame, as opposed to intrinsic, observer-indifferent features of the flow of fluid particles. For studies seeking to be consistent with observed flow physics, the latter features are relevant (Haller Reference Haller2021). This is because flow visualization experiments with dye particles highlight material (and hence objectively defined) transport barriers, which generally differ substantially from streamsurfaces in unsteady flows. At the same time, classic studies focused on momentum, energy or vorticity transport are inherently frame-dependent by the frame-dependence of all these Eulerian fields.
To reconcile these two objectives, Haller et al. (Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020) developed a theory of objectively defined barriers to the transport of active vector fields, such as the vorticity and momentum. Aksamit & Haller (Reference Aksamit and Haller2022) used this theory to locate instantaneous (Eulerian) frame-indifferent momentum transport barriers in 3-D turbulent channel flows. Active barriers turn out to be distinguished streamsurfaces of appropriately defined steady, 3-D incompressible vector fields (Haller et al. Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020). Specifically, momentum transport barriers of a Navier–Stokes velocity field 
$\boldsymbol {v}(\boldsymbol {x},t)$ at time 
$t$ are streamsurfaces of the barrier equation
\begin{equation} \boldsymbol{x}^{\prime}(s)=\boldsymbol{\Delta} \boldsymbol{v}\left(\boldsymbol{x}(s),t\right), \end{equation}
with 
$s\in {\mathbb {R}}$ denoting a parametrization of streamlines forming the streamsurfaces and 
$\boldsymbol{\Delta}$ the Laplace operator. For an incompressible velocity field, we further note that 
$\boldsymbol{\Delta} \boldsymbol {v} = -\boldsymbol {\nabla } \times \boldsymbol {\omega }$ holds.
Haller et al. (Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020) and Aksamit & Haller (Reference Aksamit and Haller2022) detected distinguished streamsurfaces of the momentum barrier equation (4.1) using active versions of some of the passive hyperbolic and elliptic LCS diagnostics reviewed by Haller (Reference Haller2015). These Lagrangian calculations involve arrays of trajectories and return diagnostic scalar fields for visual inspection rather than explicit streamsurfaces families. In the following, we will use our approach for finding approximate first integrals to obtain vortical momentum barriers as level surfaces of an approximate first integral for (4.1).
4.1. Numerical dataset
We study the 3-D incompressible, turbulent channel flow which can be found here. The friction Reynolds number is 
${Re}_{\tau } = u_{\tau }h/ \nu = 150$, where 
$u_{\tau }$ is the friction velocity, 
$h$ is the channel half-height and 
$\nu$ is the kinematic viscosity. We denote by 
$x$, 
$z$ and 
$y$ the streamwise, spanwise and wall-normal directions, respectively. The computational domain is 
$L_x = 2.5{\rm \pi} h$ long and 
$L_z = {\rm \pi}h$ wide.
The number of Fourier modes used in the simulation was 
$192$ in both the streamwise and the spanwise direction. The number of points in the wall-normal direction is 
$194$, inhomogeneously spaced so that the grid becomes more refined closer to the walls. No-slip boundary conditions were applied at the two channel walls and the governing equations were integrated forward in time with a constant pressure gradient 
$\textrm {d}p/{\textrm {d}\kern0.06em x} = -1$ enforced so as to drive the flow through the channel. Ten thousand velocity field snapshots were stored at multiples of the simulation time step 
$\boldsymbol{\Delta} t = 0.001$ once the flow reached a statistically stationary state. We will identify the 
$5000$th snapshot of this time series with the time 
$t = 0$.
4.2. Momentum barrier extraction
We seek to uncover objectively defined instantaneous momentum transport barriers corresponding to the time 
$t = 0$ as specific streamsurfaces of the vector field 
$\boldsymbol{\Delta} \boldsymbol {v}(\boldsymbol {x},0)$ in (4.1). First, as illustration of prior results by Haller et al. (Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020) on this problem, we show the FTLE field computed for 
$\boldsymbol{\Delta} \boldsymbol {v}(\boldsymbol {x},0)$ (called active FTLE or aFTLE) from a grid of 
$800 \times 1000$ initial conditions uniformly placed in the wall-normal and spanwise directions, respectively, up to 
$s = 10^{-3}$ (see figure 16a). This plot indicates the signatures of two larger barriers to momentum transport, with the first located approximately around 
$[z/h,y/h] = [0.75,1.25] \times [1.5,1.75]$. The second, a mushroom-type barrier, is located around 
$[2.25,2.75] \times [1.25,1.75]$. A plethora of other smaller-scale structures are also present in the aFTLE plots.
(a) Active FTLE (aFTLE) for the momentum barrier field (4.1) at 
$t = 0$ over a 2-D cross-section of the channel at 
$x/h = 2$. (b–d) Projections of the approximate first integral (black lines) on the same cross-section using as computational domains for the analysis the (red) boxes 
$[1,3] \times [0,2] \times [0,{\rm \pi} ]$, 
$[1,3] \times [0,2] \times [0.5,2.5]$ and 
$[1,3] \times [1.25,2] \times [0.4,1.4]$, respectively, superimposed on the aFTLE landscape.
We use three different computational domains to illustrate the potential of the algorithm presented in § 2: the boxes 
$[1,3] \times [0,2] \times [0,{\rm \pi} ]$, 
$[1,3] \times [0,2] \times [0.5,2.5]$ and 
$[1,3] \times [1.25,2] \times [0.4,1.4]$ with 
$80 \times 85 \times 80$ points in the streamwise, wall-normal and spanwise direction, respectively. In each of these boxes, the points are evenly spaced and a tri-linear interpolation scheme is employed to obtain the necessary 
$\boldsymbol{\Delta} \boldsymbol {v}$ values. Moreover, we perform another rescaling of the dummy time 
$s$ to pointwise normalize the right-hand side of (4.1). This eliminates the high norms of 
$\boldsymbol{\Delta} \boldsymbol {v}$ near the wall in the expression (2.6). If this normalization is not used, our algorithm tends to miss the turbulent regions close to the wall and touches more on the central part of the channel.
For 
$N = 13$ 
$(9170$ Fourier modes
$)$, 
$20$ isovalues of the approximate first integral are depicted in figure 16(b–d) on the 
$x/h = 2$ plane. We note the progressive improvement in the reconstruction of the first structure discussed above, as the computational domain starts to close in on it. This observation prompts us to use a greedy algorithm compartmentalizing the domain into smaller, overlapping domains. The size of these domains will be dictated by the scale of the structures to be extracted.
Our algorithm so far generates families of streamsurfaces. We now seek to locate outermost members of nested elliptical barrier families. To this end, for each of the overlapping computational boxes, we use 
$n_{iso}$ different values to produce isosurfaces in the interval 
$[ | H |_{min},| H |_{max}]$. We then group the resulting isosurfaces into foliations (see figure 20(f) for two such foliations) or, more generally, families by finding the ones that either lie entirely inside others or have an intersection volume above a certain threshold, respectively. Subsequently, we discard the foliations or families that have a number of members less than a fixed percentage of 
$n_{iso}$. We note that for the isosurface extraction, we use MATLAB's built-in function which is based on the ‘Marching cubes’ algorithm (Lorensen & Cline Reference Lorensen and Cline1987). Other techniques generating isosurfaces (see, e.g. the concept of contour trees (Carr, Snoeyink & Axen Reference Carr, Snoeyink and Axen2003)) are known to produce robust structures in a more efficient fashion and could certainly be used as viable alternatives.
We also follow Haller et al. (Reference Haller, Hadjighasem, Farazmand and Huhn2016) and define the convexity deficiency of a closed surface as the ratio of the volume between the surface and its convex hull to the volume enclosed by the surface. Here, however, we will only discard the most non-convex surfaces.
We summarize the main steps of the computations we described previously in Algorithm 1. In the next section, we will use this algorithm to uncover momentum transport barriers tied to different scales in the turbulent channel flow.
Extraction of instantaneous barriers to momentum transport
4.3. Results
4.3.1. Channel partition into large subdomains
To illustrate our algorithm, we use the entire computational domain of the direct numerical simulation (DNS) described in § 4.1. First we partition it in the following overlapping subdomains 
$[x_{1_i},x_{2_i}] \times [y_{1_j},y_{2_j}] \times [z_{1_k},z_{2_k}]$ with 
$x_1 \in \{-0.5,0.5,\ldots,5.5\}$, 
$x_2 \in \{1.5,2.5,\ldots,7.5\}$, 
$y_1 \in \{0,0.5,1\}$, 
$y_2 \in \{1,1.5,2\}$, 
$z_1 \in \{0,1,2\}$, 
$z_2 \in \{1.5,2.5,3.5\}$, 
$i = 1, \ldots, 7$, 
$j = 1, 2, 3$, 
$k = 1, 2, 3$. For each of these domains, we use a grid of 
$70 \times 75 \times 70$ evenly spaced points on which we compute the normalized velocity Laplacian. For the approximate first integral, we use 
$N=13$ (or 
$9170$ Fourier modes). This comes as a necessary trade-off between computational time and surface accuracy after we observed that the reconstructed structures do not vary significantly when we use 
$N=12$ or 
$N=14$. With it, we extract 
$n_{iso} = 40$ isosurfaces in every domain. Of them, we discard the surfaces described in steps 
$5$ and 
$6$ of Algorithm by employing 
$p_b = 2\,\%$, 
$p_f = 10\,\%$ and 
$d_{max} = 20\,\%$.
The results of this computation are shown in figure 17. We observe that the entire channel is populated by a host of different momentum transport barriers, the majority of which have a clear quasi-streamwise direction. These structures are reminiscent of those investigated in studies of wall-bounded flows (Robinson Reference Robinson1991). The structures we obtain, however, are more evenly scattered throughout the channel, appearing not only near the channel walls but also in the more quiescent central part. This is in stark contrast to typical predictions from classic vortex criteria for these flows (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988; Chong, Perry & Cantwell Reference Chong, Perry and Cantwell1990; Jeong & Hussain Reference Jeong and Hussain1995). This phenomenon, i.e. structures penetrating into and spanning the bulk flow region, has already been noted in the literature (Haller et al. Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020; Aksamit & Haller Reference Aksamit and Haller2022) by analysing the 2-D signatures of these structures via diagnostic fields. Here, however, we explicitly construct these structures via the streamsurfaces surrounding them.
Instantaneous barriers to momentum transport in the 3-D turbulent channel flow, derived using the larger partition described in § 4.3.1.
4.3.2. Channel partition into small subdomains
We now employ smaller partitions of the channel, using the overlapping subdomains 
$[x_{1_i},x_{2_i}] \times [y_{1_j},y_{2_j}] \times [z_{1_k},z_{2_k}]$ with 
$x_1 \in \{-0.5,0,\ldots,6.5\}$, 
$x_2 \in \{0.5,1,\ldots,7.5\}$, 
$y_1 \in \{0,0.5,1\}$, 
$y_2 \in \{1,1.5,2\}$, 
$z_1 \in \{0,1,2\}$, 
$z_2 \in \{1.5,2.5,3.5\}$, 
$i = 1, \ldots, 15$, 
$j = 1, 2, 3$, 
$k = 1, 2, 3$, while keeping all the other parameters the same as in § 4.3.1.
The results for this computation are shown in figure 18. In this case, we again observe momentum transport barriers that tend to align with the streamwise direction but in larger numbers, when compared with figure 17. Another notable difference is that the reconstructed structures show generally smaller scales than before. This demonstrates that our domain partition algorithm acts as a filter of various scales. This provides a natural way out of the occlusion quandary that besets many visualization approaches while retaining the algorithm's capacity to capture the smallest structures, delineated in the aFTLE plots, by using finer partitions. In the next subsection, we will demonstrate this capacity by pinpointing structures that the active diagnostics-based methods in Haller et al. (Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020) and Aksamit & Haller (Reference Aksamit and Haller2022) would likely miss. We conclude this section with figure 19 depicting a well-formed spectral gap between the smallest and the second smallest eigenvalue of 
$\boldsymbol{\mathsf{A}}$ for all the computational boxes. Finally, we carry out an unsteady barrier analysis for times varying over the interval 
$[t_{0},t_{1}] = [0,0.5]$, as presented in the supplementary movie 2. In it, the fact that some barrier surfaces seem to appear, then disappear and perhaps reappear at a later time is a direct consequence of the flat convexity deficiency threshold we allow for these surfaces (20 %) throughout the channel. We expect this phenomenon to be resolved when the purely Lagrangian approach of Haller et al. (Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020) is combined with our algorithm in follow-up work.
Instantaneous barriers to momentum transport in the 3-D turbulent channel flow, derived using the smaller partition described in § 4.3.2. The supplementary movie 2 of the supplementary material portrays the extracted structures for the time interval 
$[t_{0},t_{1}] = [0,0.5]$.
Distribution of the smallest (blue) and the second smallest (red) eigenvalues of 
$\boldsymbol{\mathsf{A}}$ over the different computational domains used for the barriers surfaces presented in figure 18.
4.3.3. Uncovering hidden structures
We now focus on the region we highlighted in figure 16(d). The signature of a vortical structure is evident from the aFTLE landscape and it is corroborated by the isolines of an approximate first integral. Zooming in on the vicinity of this structure in figure 18, however, reveals the existence of another structure (see figure 20), showing no imprint on the aFTLE landscape.
Two branches of a mushroom-like, objective vortex (in grey) captured by Algorithm 1 in the 3-D turbulent channel flow. The branches are superimposed on 2-D cross-sections of the aFTLE field at 
$x/h = 1, 1.4, 1.6, 1.9$ and 
$2.1$ in panels (a), (b), (c), (d) and (e), respectively. (f) A transverse cut at 
$x/h = 1.7$ revealing the foliations of structures constituting the two branches of the mushroom-like vortex.
To investigate this, we populate the neighbourhood around these two structures with 2-D cross-sections and compute the aFTLE field for each of them. At 
$x/h = 1$ (figure 20a), we note two mushroom-like pairs of structures, one larger and one smaller, but neither of the two display prominent closed regions that could signal a vortex. As we move towards larger 
$x/h$ values, however, the aFTLE topography changes drastically. At 
$x/h = 1.4$ (figure 20b), the larger structure starts to fall apart before it completely disintegrates in the following cross-sections. In contrast, the smaller structure develops a salient closed loop, delineated by an aFTLE ridge, which is coincident with the largest of the reconstructed surfaces. In the next two cross-sections (figure 20c,d), the smaller of the two reconstructed surfaces comes into play as we further observe both parts of the vortex pair following closely the aFTLE ridges downstream. Finally, in the last cross-section (figure 20e), the vortex pair pattern is broken up in agreement with the reconstructed vortex pair coming to a halt.
This example epitomizes two notable features of Algorithm 1. First, using finer partitions results in the reconstruction of smaller structures consistent with the well-known hierarchy of coherent structures in turbulence. Second, this reconstruction takes place in an almost automated fashion with minimal reliance on user-defined parameters and with no need to advect trajectories of the barrier equation.
5. Conclusions
We have introduced a minimization principle to find an instantaneous approximate first integral of a given 3-D vector field. Level surfaces of this approximate first integral are expected to provide automated approximations to vortical streamsurfaces highlighting the instantaneous elliptic regions of the vector field. We have also proposed and tested a numerical algorithm for solving this minimization problem by means of total least squares and ridge regression. This algorithm has performed well on various solutions of the steady 3-D Euler equations on triply periodic and non-periodic domains. For steady velocity fields, the parametrized family of streamsurfaces rendered by our approach approximate elliptic (toroidal or cylindrical) LCSs without the need to advect a large number of trajectories required by Lagrangian methods (Haller et al. Reference Haller, Hadjighasem, Farazmand and Huhn2016). We have also illustrated on a 3-D junction flow the applicability of our algorithm for velocity fields defined as numerical datasets produced from computational fluid dynamics simulations.
We have additionally used the described algorithm to extract objectively defined instantaneous momentum transport barriers in a 3-D turbulent channel flow in an almost-automated fashion. Such barriers are streamsurfaces of an incompressible barrier equation defined by the Laplacian of the velocity field (Haller et al. Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020; Aksamit & Haller Reference Aksamit and Haller2022). We have extracted vortical momentum transport barriers without the need to advect arrays of trajectories and found such barriers across multiple spatial scales.
The generality of the presented algorithm provides a fertile ground for its use in the visualization of the important features in a number of diverse applications. A straightforward example is the barrier equation derived by Haller et al. (Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020) for the transport of vorticity. Furthermore, the purely Lagrangian approach of Haller et al. (Reference Haller, Katsanoulis, Holzner, Frohnapfel and Gatti2020) could be combined with our algorithm to produce LCSs for unsteady simulations. This method, unlike the instantaneous approach used for the production of supplementary movie 2, will yield smoothly varying structures that are also experimentally verifiable. Another application of our algorithm could be the visualization of the magnetic field resulting from magnetohydrodynamic turbulence simulations (Biskamp Reference Biskamp2003). Future research can also apply the algorithm to direction fields like those emerging from studies of 3-D LCSs (Oettinger & Haller Reference Oettinger and Haller2016), to 3-D material barriers to diffusive transport (Haller, Karrasch & Kogelbauer Reference Haller, Karrasch and Kogelbauer2018) or to uncover LCSs acting as barriers in diverse applications (Martínez et al. Reference Martínez2021). Such studies could extend the automated extraction of 2-D LCSs and transport barriers (Katsanoulis & Haller Reference Katsanoulis and Haller2019; Katsanoulis et al. Reference Katsanoulis, Farazmand, Serra and Haller2020) to 3-D flows, facilitating their integration to numerical simulation codes.
Alternative solutions to the same optimization problem for an approximate first integral are certainly viable to investigate in future work. These could involve different basis functions (such as Chebyshev polynomials, wavelets or radial basis functions) instead of Fourier expansions. Techniques to mitigate the Gibbs phenomenon tied to non-periodic domains (Gegenbauer polynomials) (Gottlieb & Shu Reference Gottlieb and Shu1997) are also feasible to consider. Similarly, the effect of using Chebyshev nodes, instead of the uniform grids we considered here, to reduce the impact Runge's phenomenon has near the boundaries could be examined.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2022.992.
Funding
S.K. and G.H. acknowledge financial support from the Priority Program SPP 1881 (Turbulent Superstructures) of the German National Science Foundation (DFG). We are grateful to Dr D. Gatti for providing us with the 3-D turbulent channel flow dataset he originally generated as a benchmark case for the same program.
Declaration of interests
The authors report no conflict of interest.
Code and data availability
Example scripts generating the approximate first integral for the ABC flow are available at https://github.com/katsanoulis/Approximate_First_Integral/. The V-junction and turbulent channel flow datasets as well as their respective surface extraction scripts are available upon request from the first author.
Appendix A. Typical flow visualization techniques
In this section, we apply classic vortex visualization methods in some of the flows for which we have constructed approximate first integrals. First, we test two isosurface-based criteria that are widely used in the literature to locate vortical structures, i.e. the 
$\lambda _2$-criterion of Jeong & Hussain (Reference Jeong and Hussain1995) and the 
$Q$-criterion of Hunt et al. (Reference Hunt, Wray and Moin1988). Both of them are defined from the decomposition
\begin{equation} \boldsymbol{\nabla} \boldsymbol{v} = \boldsymbol{\mathsf{S}} + {\boldsymbol{\mathsf{\Omega}}}, \end{equation}
where 
$\boldsymbol{\mathsf{S}}=[\boldsymbol {\nabla } \boldsymbol {v} + (\boldsymbol {\nabla } \boldsymbol {v})^\top ]/2$ is the rate-of-strain tensor and 
${\boldsymbol{\mathsf{\Omega}}}=[\boldsymbol {\nabla } \boldsymbol {v} - (\boldsymbol {\nabla } \boldsymbol {v})^\top ]/2$ is the vorticity tensor.
According to the 
$\lambda _2$-criterion, vortical regions coincide with domains where
\begin{equation} \lambda_2(\boldsymbol{\mathsf{S}}^2 + {\boldsymbol{\mathsf{\Omega}}}^2) < 0 \end{equation}
with 
$\lambda _2(\boldsymbol{\mathsf{B}})$ denoting the intermediate eigenvalue of the symmetric tensor 
$\boldsymbol{\mathsf{B}}$. Similarly, the 
$Q$-criterion identifies vortical regions as those where
\begin{equation} Q = \frac{1}{2} \left[ \lVert {\boldsymbol{\mathsf{\Omega}}} \rVert^2 - \lVert \boldsymbol{\mathsf{S}} \rVert^2 \right] > 0 \end{equation}
with 
$\lVert \boldsymbol{\mathsf{B}} \rVert$ denoting the Frobenius norm of 
$\boldsymbol{\mathsf{B}}$.
We stress here that both of these criteria are only Galilean-invariant and, thus, depend on the frame of reference they are used on. Of equal importance is that, according to their definitions, they highlight vortical regions rather than surfaces. To bypass this, a specific value is usually chosen and the resulting isosurfaces are identified as the so-called vortices. Such a choice, however, would be justifiable, if the chosen isovalue was close to zero, as the original definitions postulated. Unfortunately, in typical flow visualizations, these isovalues are tuned to arbitrary values to match user expectations regarding the vortical features.
We highlight this in figure 21 using the 
$\lambda _2$-criterion for the non-integrable ABC flow. Specifically, we use one value close to zero (figure 21a) and one which corresponds to a drastically different topology for the resulting structures (figure 21b). We note that both of these values fail to capture even a single family of invariant tori. Indeed, the reconstructed surfaces tend to either have holes in the vicinity of the invariant tori (figure 21a) or be misaligned with the vortex cores (figure 21b). This is in stark contrast to the approximate-first-integral-based tori we have found in § 3.1.2.
$\lambda _2$-criterion-based isosurfaces against Poincaré maps in the non-integrable ABC flow. The isosurfaces correspond to (a) 
$\lambda _2=-0.02$ and (b) 
$\lambda _2=-2.4$. The Poincaré maps are computed from a grid of 
$20 \times 20$ initial conditions on the 
$x=0$, 
$y=0$ and 
$z=0$ planes running up to arclength 
$10^4$.
A similar conclusion can be drawn for the flow inside the V-junction that we tested in § 3.4. Our approximate-first-integral algorithm correctly pinpoints the exact position of the recirculation bubbles in this flow. In contrast, we show in figure 22 the drastically different vortical features we obtain for different isovalues of the 
$Q$-criterion. Despite the four recirculation bubbles that are formed downstream, the isosurfaces based on the 
$Q$-criterion locate invariably two vortical features whose length also varies substantially for different 
$Q$ values. What is even more remarkable is that such 
$Q$-criterion-based vortical features persist even for similar flows with slightly different Reynolds numbers that exhibit no fluid recirculation at all (see figure 22c). Importantly, our method finds no structure in the same spatial domain for these Reynolds numbers and hence correctly signals the lack of a recirculation bubble.
(a,b) 
$Q$-criterion-based isosurfaces against the recirculation bubbles for the flow inside the V-junction of § 3.4 
$({Re} = 230)$. The isosurfaces correspond to (a) 
$Q=0.02$ and (b) 
$Q=50$. (c) 
$Q$-criterion-based isosurfaces for 
$Q = 50$ for a perturbed solution 
$({Re} = 180)$ where no recirculation bubbles are formed.
Finally, in figure 23(a,b), we present representative streamlines for the turbulent channel flow of § 4 based on the barrier field 
$\boldsymbol{\Delta} \boldsymbol {v} (\boldsymbol {x})$, the velocity field 
$\boldsymbol {v} (\boldsymbol {x})$ in figure 23(c,d) and the vorticity field 
$\boldsymbol \omega (\boldsymbol {x})$ in figure 23(e,f) at time 
$t=0$. For all three fields, we choose the same initial conditions, i.e. points that approximately lie on the ridge of the aFTLE field (computed based on the 
$\boldsymbol{\Delta} \boldsymbol {v} (\boldsymbol {x})$ field) we identified in figure 20(c) at 
$x/h=1.6$.
Streamlines of different vector fields emanating from points lying on the ridge of the aFTLE field at 
$x/h = 1.6$. Different views of streamlines of the (a,b) momentum barrier field, (c,d) velocity field and (e,f) vorticity field at 
$t=0$.
Based on these calculations, we draw the following conclusions. First, we observe the barrier streamlines being wrapped around the approximate-first-integral-based structure following even small protrusions on its surface like the one seen in figure 23(a). This structure is, thus, a correct approximation to an invariant manifold of the barrier field in agreement with the best-fit streamsurface of figure 24(a,b) as well as the aFTLE landscape of figure 20. For longer integration dummy times (figure 24c,d), however, we observe that this delineation of the barrier surface using streamlines quickly comes to a halt. Indeed, the streamlines fail to capture a significant portion of the barrier surface as 
$x/h$ becomes larger, whereas for smaller 
$x/h$, they develop convoluted patterns that eventually fall apart resulting in elongated streamsurfaces that show no imprint on the aFTLE landscape. This constitutes the main reason why techniques aiming at the reconstruction of exact streamsurfaces through a better seed placement are not well suited for such flows.
Different views of streamlines of the momentum barrier field originating from points lying on the ridge of the aFTLE field presented in figure 23. The integration in panels (a,b) and (c,d) corresponds to different dummy final times 
$s_{1,(a,b)}$ and 
$s_{1,(c,d)}$ with 
$s_{1,(a,b)} < s_{1,(c,d)}$. The best-fit streamsurfaces are depicted in red and are compared against the approximate-first-integral-based structure in grey.
Second, we note that the velocity streamlines form a misaligned (with respect to the extracted structure) tube that gives no indication of a vortical feature. Third, we remark that the vorticity streamlines correctly outline the outer shape of the mushroom-shaped structure imprinted on the aFTLE landscape without a detailed delineation, however, of the two branches that constitute it.
Appendix B. Least squares for homogeneous systems
Theorem B.1 Let 
$\boldsymbol{\mathsf{C}}=\boldsymbol{\mathsf{U}}{\boldsymbol{\mathsf{\Sigma}}} \boldsymbol{\mathsf{V}}^{*}$ be the singular value decomposition of a matrix 
$\boldsymbol{\mathsf{C}}$. Let, also, 
$\boldsymbol {v}_{1}, \ldots, \boldsymbol {v}_k$ be the last 
$k$ columns of 
$\boldsymbol{\mathsf{V}}$ whose corresponding singular values are equal to the smallest singular value 
$\sigma _1$. Then, all the linear combinations of the form
\begin{equation} \boldsymbol{x} = c_1 \boldsymbol{v}_{1} + \cdots + c_k \boldsymbol{v}_k \end{equation}with
\begin{equation} c_1^2 + \cdots + c_k^2 = 1 \end{equation}are unit-norm least-squares solutions to the homogeneous system
\begin{equation} \boldsymbol{\mathsf{C}}\boldsymbol{x} = \boldsymbol{0}. \end{equation}Proof. We want to find the solution 
$\boldsymbol {x}$ with 
$|\boldsymbol {x} | = 1$ that minimizes 
$| \boldsymbol{\mathsf{C}} \boldsymbol {x} |$ or 
$| \boldsymbol{\mathsf{U}} {\boldsymbol{\mathsf{\Sigma}}} \boldsymbol{\mathsf{V}}^{*} \boldsymbol {x}|$. Because 
$\boldsymbol{\mathsf{U}}$ is unitary and, thus, acts as an isometry, this is equivalent to minimizing 
$| {\boldsymbol{\mathsf{\Sigma}}} \boldsymbol{\mathsf{V}}^{*} \boldsymbol {x}|$ or 
$|{\boldsymbol{\mathsf{\Sigma}}} \boldsymbol {y}|$ with 
$\boldsymbol {y} = \boldsymbol{\mathsf{V}}^{*} \boldsymbol {x}$. Here, 
$\boldsymbol{\mathsf{V}}$ is also unitary, so 
$|\boldsymbol {y} | = 1$ is equivalent to 
$|\boldsymbol {x} | = 1$. Therefore we want the unit-norm vector 
$\boldsymbol {y}$ that minimizes 
$|{\boldsymbol{\mathsf{\Sigma}}} \boldsymbol {y}|$, i.e. the quantity
\begin{equation} \sigma_1^2 y_1^2 + \cdots + \sigma_n^2 y_n^2. \end{equation}By the definition of SVD, we have
\begin{equation} \left.\begin{array}{c@{}} \sigma_{1}\geq\sigma_{1}\Leftrightarrow\sigma_{1}^{2}\geq\sigma_{1}^{2}\\ \vdots\\ \sigma_{n}\geq\sigma_{1}\Leftrightarrow\sigma_{n}^{2}\geq\sigma_{1}^{2} \end{array}\right\} \left.\Leftrightarrow\begin{array}{c@{}} \sigma_{1}^{2}y_{1}^{2}\geq\sigma_{1}^{2}y_{1}^{2}\\ \vdots\\ \sigma_{n}^{2}y_{n}^{2}\geq\sigma_{1}^{2}y_{n}^{2} \end{array}\right\} \end{equation}and after summing up all the inequalities, we obtain
\begin{equation} \sigma_{1}^{2}y_{1}^{2}+\sigma_{2}^{2}y_{2}^{2}+\cdots+\sigma_{n}^{2}y_{n}^{2}\geq\left(y_{1}^{2}+y_{2}^{2}+\cdots+y_{n}^{2}\right)\sigma_{1}^{2}=\sigma_{1}^{2}, \end{equation}with the equality holding if
\begin{equation} y_{n-k+1}=y_{n-k+2}=\cdots=y_{n}=0. \end{equation}
We then have 
$\boldsymbol {x} = \boldsymbol{\mathsf{V}} \boldsymbol {y} = y_1 \boldsymbol {v}_1 + \cdots + y_n \boldsymbol {v}_n$ and, thus, (B7) is equivalent to
\begin{equation} \boldsymbol{x} = c_1 \boldsymbol{v}_{1} + \cdots + c_k \boldsymbol{v}_k \end{equation}
with 
$c_1 = y_{1}, \dotsc, c_k=y_k$ and 
$c_1^2 + \cdots + c_k^2 = 1$.
Appendix C. Solving for H after imposing the real constraint
We define a complex vector 
$\boldsymbol {h} \in \mathbb {C}^n$, where 
$n$ represents the number of vectors 
$\boldsymbol {k} \in \mathbb {Z}^3$ such that 
$|\boldsymbol {k}|\leq N$ and we also represent the individual entries of 
$\boldsymbol {h}$ as 
$h_{\boldsymbol {k}}$. Further, in conjunction with § 2, we define a 
$m\times n$ matrix, where 
$m$ represents the total number of grid points. We are interested in minimizing 
$\boldsymbol {h}^{*} \boldsymbol{\mathsf{C}}^{*} \boldsymbol{\mathsf{C}} \boldsymbol {h}$ subject to the constraints 
$\boldsymbol {h}^{*} \boldsymbol {h} = 1$ and 
$h_{-\boldsymbol {k}}=h^*_{\boldsymbol {k}}$. The second constraint is not natural to solve using known optimization techniques, due to which we will modify the above problem rearranging 
$\boldsymbol {h}$ as a vector in 
$\mathbb {R}^{2n}$.
Naively, one can construct a column vector with first 
$n$ real entries of 
$\boldsymbol {h}$ followed by their corresponding imaginary entries
\begin{equation} \boldsymbol{v} = \begin{pmatrix} \operatorname{Re}(\boldsymbol{h}) \\ \operatorname{Im}(\boldsymbol{h}) \end{pmatrix}. \end{equation}
Similarly, one can reorganize 
$\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{C}}^{*} \boldsymbol{\mathsf{C}}$ and obtain
\begin{equation} \boldsymbol{\mathsf{D}} = \begin{bmatrix} \boldsymbol{\mathsf{A}}_r & -\boldsymbol{\mathsf{A}}_{\rm i} \\ \boldsymbol{\mathsf{A}}_{\rm i} & \boldsymbol{\mathsf{A}}_r \end{bmatrix}, \end{equation}
where 
$\boldsymbol{\mathsf{A}}_r = ({\boldsymbol{\mathsf{A}}+\boldsymbol{\mathsf{A}}^\textrm {T}})/{2}$ and 
$\boldsymbol{\mathsf{A}}_{\rm i} = ({\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{A}}^\textrm {T}})/{2{\rm i}}$. Thus, we have 
$\boldsymbol {v}^\textrm {T}\boldsymbol{\mathsf{D}}\boldsymbol {v}=\boldsymbol {h}^{*}\boldsymbol{\mathsf{A}}\boldsymbol {h}$ subject to the constraints 
$\boldsymbol {v}^\textrm {T}\boldsymbol {v}=\boldsymbol {h}^{*} \boldsymbol {h} = 1$ and 
$\boldsymbol{\mathsf{U}}\boldsymbol {v}=0$. In our case, 
$\boldsymbol{\mathsf{U}}$ is a matrix of the form
\begin{gather} \boldsymbol{\mathsf{U}} = \begin{bmatrix} \boldsymbol{\mathsf{K}}_{(n\times n)} & \boldsymbol{\mathsf{J}}_{(n\times n)} \end{bmatrix}_{n \times 2n} \end{gather}
\begin{gather}{\mathsf{K}}_{ij} = \begin{cases} 1 & i=j,i< n/2,j< n/2 \\ -1 & i=n/2+j, i< n/2,j< n/2\\ 0 \end{cases} \bigg| {\mathsf{J}}_{ij} = \begin{cases} 1 & i=j,i< n/2,j< n/2 \\ 1 & i=n/2+j, i< n/2,j< n/2\\ 0 \end{cases} \end{gather}The above optimization problem has a solution by Golub (Reference Golub1973). We state a form of the theorem that will be useful for us.
Theorem C.1 Minimizing 
$\boldsymbol {v}^\textrm {T} \boldsymbol{\mathsf{D}}\boldsymbol {v}$ subject to 
$\boldsymbol {v}^\textrm {T}\boldsymbol {v}=1$, 
$\boldsymbol{\mathsf{U}}\boldsymbol {v}=0$, where 
$\boldsymbol {v}\in \mathbb {R}^{2n}$ is equivalent to minimizing 
$\boldsymbol {l}^\textrm {T} \boldsymbol{\mathsf{E}} \boldsymbol {l}$ subject to 
$\boldsymbol {l}^\textrm {T}\boldsymbol {l}=1$, where 
$\boldsymbol {l}\in \mathbb {R}^n$ provided 
$\boldsymbol{\mathsf{U}}$ is a 
$n\times 2n$ matrix.
Proof. Following Golub (Reference Golub1973), we can write
\begin{equation} \boldsymbol{\mathsf{U}}^{\rm T} = \boldsymbol{\mathsf{QR}} = \boldsymbol{\mathsf{Q}}_{2n\times 2n}\begin{bmatrix} \boldsymbol{\mathsf{R}}'_{n\times n} \\ \boldsymbol{0}_{n\times n} \end{bmatrix}_{2n\times n}, \end{equation}
where 
$\boldsymbol{\mathsf{Q}}$ is an orthogonal matrix and 
$\boldsymbol{\mathsf{R}}'$ is an upper triangular matrix. The constraint then reads
\begin{gather} \boldsymbol{\mathsf{R}}^{\rm T}\boldsymbol{\mathsf{Q}}^{\rm T}\boldsymbol{v}=0, \end{gather}
\begin{gather}\boldsymbol{\mathsf{R}}^{\rm T}\boldsymbol{v}_{\boldsymbol{Q}}=0. \end{gather}
The vector 
$\boldsymbol {v}_{\boldsymbol{Q}} = \left [\begin {smallmatrix} \boldsymbol {0}_{n\times 1} \\ \boldsymbol {l}_{n \times 1}\end {smallmatrix}\right ]$ is a solution to the above constraint. Expressing the objective function in terms of 
$\boldsymbol {v}_{\boldsymbol{Q}}$, we have
\begin{equation} \boldsymbol{v}^{\rm T}\boldsymbol{\mathsf{D}}\boldsymbol{v} = \boldsymbol{v}^{\rm T}\boldsymbol{\mathsf{QQ}}^{{-}1}\boldsymbol{\mathsf{D}}(\boldsymbol{\mathsf{Q}}^{{-}1})^{\rm T}\boldsymbol{\mathsf{Q}}^{\rm T}\boldsymbol{v}=\boldsymbol{v}^{\rm T}_{\boldsymbol{Q}} \begin{bmatrix}(\boldsymbol{\mathsf{Q}}^{\rm T}\boldsymbol{\mathsf{DQ}})_{11} & (\boldsymbol{\mathsf{Q}}^{\rm T}\boldsymbol{\mathsf{DQ}})_{12} \\ (\boldsymbol{\mathsf{Q}}^{\rm T}\boldsymbol{\mathsf{DQ}})_{21} & (\boldsymbol{\mathsf{Q}}^{\rm T}\boldsymbol{\mathsf{DQ}})_{22} \end{bmatrix}\boldsymbol{v}_{\boldsymbol{Q}} = \boldsymbol{l}^{\rm T}\boldsymbol{\mathsf{E}}\boldsymbol{l} \end{equation}
subject to the constraint 
$\boldsymbol {l}^\textrm {T}\boldsymbol {l}=1$ with 
$\boldsymbol{\mathsf{E}} = (\boldsymbol{\mathsf{Q}}^\textrm {T}\boldsymbol{\mathsf{DQ}})_{22}$.
Using this theorem, we can see that the eigenvector of 
$\boldsymbol{\mathsf{E}}$ corresponding to the minimum eigenvalue would be the solution (since 
$\boldsymbol{\mathsf{E}}$ is a symmetric matrix). The vector of interest 
$\boldsymbol {h}$ can be reconstructed as follows:
\begin{equation} \boldsymbol{h} = [\boldsymbol{\mathsf{Q}}(\boldsymbol{v}_{\boldsymbol{Q}})_{min}]_{1} + {\rm i} [\boldsymbol{\mathsf{Q}}(\boldsymbol{v}_{\boldsymbol{Q}})_{min}]_2,\boldsymbol{\mathsf{Q}}(\boldsymbol{v}_{\boldsymbol{Q}})_{min} = \begin{pmatrix} \{[\boldsymbol{\mathsf{Q}}(\boldsymbol{v}_{\boldsymbol{Q}})_{min}]_1\}_{n\times 1} \\ \{[\boldsymbol{\mathsf{Q}}(\boldsymbol{v}_{\boldsymbol{Q}})_{min}]_2\}_{n\times 1} \end{pmatrix}. \end{equation}The results of this solution method for the non-integrable ABC flow on three different planes are depicted in figure 25. Again, we note the very good agreement of the reconstructed level sets with both the Poincaré maps and the level sets of figure 3.
Analysis of the non-integrable ABC flow using a computational grid of 
$100^3$ points and 
$9170$ Fourier modes after constraining 
$H$ to be a real scalar field. Level sets of the reconstructed first integral at (a) 
$z=0$, (b) 
$y=0$ and (c) 
$x=0$. The Poincaré map is overlaid on each section based on a uniform grid of 
$20 \times 20$ initial conditions.
Appendix D. Efficient computation of SVD for tall-skinny matrices
Let 
$\boldsymbol{\mathsf{C}}$ be an 
$m \times n$ matrix with 
$m \gg n$. To make use of the tall-thin structure of this matrix, we follow an approach similar to the one described by Schmidt (Reference Schmidt2020) and factor
\begin{equation} \boldsymbol{\mathsf{C}}_{m\times n}=\boldsymbol{\mathsf{Q}}_{m\times n} \boldsymbol{\mathsf{R}}_{n\times n} \end{equation}
using the so-called thin QR factorization of 
$\boldsymbol{\mathsf{C}}$. The SVD of 
$\boldsymbol{\mathsf{R}}$ gives
\begin{equation} \boldsymbol{\mathsf{C}}=\boldsymbol{\mathsf{Q}} \boldsymbol{\mathsf{R}}=\boldsymbol{\mathsf{Q}}(\boldsymbol{\mathsf{U}}_{R} {\boldsymbol{\mathsf{\Sigma}}}_{R} \boldsymbol{\mathsf{V}}_{R}^{{\rm T}})=(\boldsymbol{\mathsf{Q}} \boldsymbol{\mathsf{U}}_{R}) {\boldsymbol{\mathsf{\Sigma}}}_{R} \boldsymbol{\mathsf{V}}_{R}^{{\rm T}}. \end{equation}
This is in turn a singular value decomposition of 
$\boldsymbol{\mathsf{C}}$ because 
$\boldsymbol{\mathsf{Q}} \boldsymbol{\mathsf{U}}_{R}$ is unitary as the product of unitary matrices. To avail ourselves of this result, we can split 
$\boldsymbol{\mathsf{C}} =\left [\begin {smallmatrix} \boldsymbol{C}_{1}\\ \boldsymbol{C}_{2}\end {smallmatrix}\right ]$, followed by 
$\boldsymbol{\mathsf{C}}_{i}=\boldsymbol{\mathsf{Q}}_{i}\boldsymbol{\mathsf{R}}_{i} \ (i=1,2)$ and, thus, we obtain
\begin{equation} \boldsymbol{\mathsf{C}}=\left[\begin{array}{@{}c@{}c@{}} \boldsymbol{\mathsf{Q}}_{1} & 0\\ 0 & \boldsymbol{\mathsf{Q}}_{2} \end{array}\right]\left[\begin{array}{c} \boldsymbol{\mathsf{R}}_{1}\\ \boldsymbol{\mathsf{R}}_{2} \end{array}\right]. \end{equation}
This is not yet a QR factorization of 
$\boldsymbol{\mathsf{C}}$ as the right factor is not upper triangular. We then perform another thin QR factorization on 
$\left [\begin {smallmatrix} \boldsymbol{\mathsf{R}}_{1}\\ \boldsymbol{\mathsf{R}}_{2} \end {smallmatrix}\right ]$ and obtain
\begin{equation} \boldsymbol{\mathsf{C}}=\left[\begin{array}{@{}c@{}c@{}} \boldsymbol{\mathsf{Q}}_{1} & 0\\ 0 & \boldsymbol{\mathsf{Q}}_{2} \end{array}\right] \boldsymbol{\mathsf{Q}} \boldsymbol{\mathsf{R}}, \end{equation}
which is a QR factorization of 
$\boldsymbol{\mathsf{C}}$. Based on that, we can partition the original computation domain into smaller sub-domains, compute the 
$\boldsymbol{\mathsf{R}}_i$ matrices for each of them, combine these to produce the final 
$\boldsymbol{\mathsf{R}}$ matrix as in (D4) and, finally, run SVD on this 
$n \times n$ matrix. This procedure is readily parallelizable.
For the purposes of our computation described in § 3.1.2, we used 
$10$ cores on the Euler cluster of ETH Zurich and within each core, we chose 
$5$ more partitions. This resulted in approximately the same memory footprint as for the solution of the eigenvalue problem. If the memory is not enough, the number of partitions can be increased with a subsequent increase in the computational time.
Appendix E. A different solution to the optimization problem
Let us assume that we start by fixing the Fourier coefficient of particular modes. To avoid inducing unnecessary inhomogeneity to our solution, we are going to set 
$\hat {H}_{\boldsymbol {k}}=1$ for 
$\boldsymbol {k}=(1,0,0)$, 
$\boldsymbol {k}=(0,1,0)$, 
$\boldsymbol {k}=(0,0,1)$. Similarly, we impose 
$\hat {H}_{\boldsymbol {-k}}=\hat {H}_{\boldsymbol {k}}^{*}=1$. We further denote by
$$\begin{gather} \mathcal{K}= \{ \boldsymbol{k} |\boldsymbol{k}\in \mathbb{Z}^{3}\land |\boldsymbol{k}| \leq N \land \boldsymbol{k}\neq(1,0,0)\land\boldsymbol{k}\neq({-}1,0,0)\land\boldsymbol{k}\neq(0,1,0)\nonumber\\ \land \boldsymbol{k}\neq(0,-1,0)\land\boldsymbol{k}\neq(0,0,1)\land\boldsymbol{k}\neq(0,0,-1) \} \end{gather}$$the set of the remaining Fourier modes. Then, the invariance condition reads
\begin{equation} \boldsymbol{\mathsf{C}}\boldsymbol{\boldsymbol{h}}=\boldsymbol{b}, \end{equation}
where 
${\mathsf{C}}_{ij}=\exp ({\textrm {i}\boldsymbol {k}_{j}\boldsymbol {\cdot}\boldsymbol {x}_{i}})\boldsymbol {k}_{j}\boldsymbol {\cdot}\boldsymbol {v}_{i}$, 
$\boldsymbol {\boldsymbol {h}}=\{ \hat {H}_{\boldsymbol {k}}|\boldsymbol {k}\in \mathbb {Z}^{3}\land \boldsymbol {k}\in \mathcal {K}\}$ and 
$\boldsymbol {b}=[b_{1} \quad b_{2} \quad \cdots \quad b_{l}$ 
$\cdots \quad b_{m}]^{\textrm {T}}$ with
\begin{equation} b_{l}={-}2{\rm i}\left(\sin x{v}_{x}+\sin y{v}_{y}+\sin z{v}_{z}\right). \end{equation}Essentially, what we have attained here is to transform the homogeneous system of equations to an inhomogeneous one for which we can now use ordinary least squares or ridge regression. The ridge-regression solution is given as
\begin{equation} \hat{\boldsymbol{h}} = \boldsymbol{\mathsf{V}} \left( {\boldsymbol{\mathsf{\Sigma}}}^{\top} {\boldsymbol{\mathsf{\Sigma}}} + \lambda \boldsymbol{\mathsf{I}}\right)^{{-}1} {\boldsymbol{\mathsf{\Sigma}}}^{\top} \boldsymbol{\mathsf{U}}^{*} \boldsymbol{b}, \end{equation}
whereas for 
$\lambda = 0$, we obtain the ordinary least-squares solution. This expression shows that the trick we employed in the homogeneous case with the QR decomposition cannot be applied here given that the solution depends on the explicit construction of 
$\boldsymbol{\mathsf{U}}$. We are, therefore, limited with regard to the maximum size of the computational grid we can use before we run into memory issues.
Nonetheless, in figure 26, we provide the two solutions for the non-integrable ABC flow using 
$60$ points per direction and the same number of Fourier modes as before. For the ridge regression, we used cross-validation and kept the solution with the smallest least-squares error out of the solutions generated for 
$\lambda \in \{1,10^{-1},10^{-2},10^{-3},10^{-4}\}$. The results appear to be inferior to the homogeneous-system-based solutions suggesting that this approach is recommended only for numerical datasets defined over smaller grids.
Analysis of the non-integrable ABC flow using a computational grid of 
$60^3$ points and 
$9170$ Fourier modes. Level sets of the reconstructed first integral at (a,d) 
$z=0$, (b,e) 
$y=0$ and (c,f) 
$x=0$. Panels (a–c) are constructed using ordinary least squares, whereas panels (d–f) are generated using ridge regression. The Poincaré map is overlaid on each section based on a uniform grid of 
$20 \times 20$ initial conditions.
Open access
