3.1 Definitions of primal and tangent solutions

Both CLVs and shadowing directions are tangent solutions, which describe evolutions of trajectory perturbations caused by perturbations on a dynamical system. There are two kinds of perturbations we can perform on a dynamical system given in (3.1): perturbation on initial conditions and on system parameters. The corresponding tangent solutions are homogeneous and inhomogeneous tangent solutions, respectively. As we shall see, CLVs are special homogeneous tangent solutions, while shadowing directions are special inhomogeneous tangent solutions.

We write the governing equation of the flow field in the form of a general dynamical system, which will be referred to as the primal system:

(3.1) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{u}}{\text{d}t}=\boldsymbol{f}(\boldsymbol{u},s),\quad \boldsymbol{u}|_{t=0}=\boldsymbol{u}^{0}+\boldsymbol{v}^{0}s+\boldsymbol{w}^{0}\unicode[STIX]{x1D719}.\end{eqnarray}$$

Here $\boldsymbol{f}(\boldsymbol{u},s):X\times \mathbb{R}\rightarrow X$ is a smooth function, $\boldsymbol{u}$ are state variables, and $s$ is the system parameter. The initial condition is $\boldsymbol{u}^{0}$ with a possible perturbation in the direction of $\boldsymbol{v}^{0}$ or $\boldsymbol{w}^{0}$ , controlled by $s$ and another parameter $\unicode[STIX]{x1D719}$ . A primal solution $\boldsymbol{u}:\mathbb{R}\rightarrow X$ maps time to a compact phase space $X$ , where a single point represents an entire 3D flow field at an instant; thus $\boldsymbol{u}(t)$ can as well be imagined as a trajectory in  $X$ .

We assume the system is autonomous, that is, $\boldsymbol{f}$ does not depend on time. We could extend our discussions to some special non-autonomous cases, such as those with some compactness in the time dependence, so that the ergodic theory is still available. For example, for cases whose time dependence can be decomposed into finitely many frequencies, we can extend the phase space by adding finitely many tori and still maintain a compact phase space. On this extended phase space our analysis will still be valid; for example, Antali & Horvath (Reference Antali and Horvath2018) used the shadowing method to compute sensitivities. However, our discussions in this paper do not extend to non-autonomous systems in the most general setting.

Our dynamical system considered in this paper is given by the semi-discretized Navier–Stokes equation without turbulence models. More specifically, on the spatial direction, the flow field, together with boundary conditions, is represented by the finite-volume method on our meshes. The state variables are conservative variables, that is, $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D70C}\boldsymbol{U}$ and $\unicode[STIX]{x1D70C}E$ . Hence our phase space $X\subset \mathbb{R}^{m}$ is finite-dimensional, with $m$ being the DOFs. We further assume that $X$ is a bounded hence compact subset. For the coarser mesh, the number of DOFs is approximately $1.9\times 10^{6}$ ; for the finer mesh, it is approximately $3.8\times 10^{6}$ . On the other hand, we regard the temporal direction as being continuous; thus our dynamical system is continuous, and later the theories we use are for continuous dynamics.

We can differentiate (3.1) with respect to $\unicode[STIX]{x1D719}$ , and define $\boldsymbol{w}=\text{d}\boldsymbol{u}/\text{d}\unicode[STIX]{x1D719}$ . Then $\boldsymbol{w}$ satisfies the so-called homogeneous tangent equation:

(3.2) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{w}}{\text{d}t}-\unicode[STIX]{x2202}_{\boldsymbol{u}}\,\boldsymbol{f}\boldsymbol{w}=0,\quad \boldsymbol{w}(t=0)=\boldsymbol{w}^{0},\end{eqnarray}$$

where $\unicode[STIX]{x2202}_{\boldsymbol{u}}\,\boldsymbol{f}\in \mathbb{R}^{m\times m}$ , and the ordinary differential equation (ODE) is called the homogeneous tangent ODE. Here $\boldsymbol{w}$ reflects the trajectory perturbation caused by perturbing the initial condition in the direction of $\boldsymbol{w}^{0}$ .

On the other hand, we can differentiate (3.1) with respect to $s$ , and define $\boldsymbol{v}=\text{d}\boldsymbol{u}/\text{d}s$ . Then $\boldsymbol{v}$ satisfies the inhomogeneous tangent equation:

(3.3) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{v}}{\text{d}t}-\unicode[STIX]{x2202}_{\boldsymbol{u}}\,\boldsymbol{f}\boldsymbol{v}=\unicode[STIX]{x2202}_{s}\,\boldsymbol{f},\quad \boldsymbol{v}(t=0)=\boldsymbol{v}^{0},\end{eqnarray}$$

where $\unicode[STIX]{x2202}_{s}\,\boldsymbol{f}\in \mathbb{R}^{m}$ , and the ODE is called the inhomogeneous tangent ODE. Here $\boldsymbol{v}$ reflects the trajectory perturbation caused by perturbing the system parameter $s$ , which not only affects the governing differential equation, but also potentially impacts the initial condition in the direction of $\boldsymbol{v}^{0}$ . The solution set of an inhomogeneous tangent ODE can be written as a particular inhomogeneous $\boldsymbol{v}^{\ast }$ plus the solution set of the homogeneous tangent ODE, which is a linear subspace due to homogeneity.

3.2 Definitions of covariant Lyapunov vectors and hyperbolicity

A CLV $\unicode[STIX]{x1D73B}(t)$ is a homogeneous tangent solution whose norm behaves like an exponential function of time. That is, there are $C_{1},C_{2}>0$ and $\unicode[STIX]{x1D706}\in \mathbb{R}$ such that, for any $t$ ,

(3.4) $$\begin{eqnarray}C_{1}\text{e}^{\unicode[STIX]{x1D706}t}\Vert \unicode[STIX]{x1D73B}(\boldsymbol{u}(0))\Vert \leqslant \Vert \unicode[STIX]{x1D73B}(\boldsymbol{u}(t))\Vert \leqslant C_{2}\text{e}^{\unicode[STIX]{x1D706}t}\Vert \unicode[STIX]{x1D73B}(\boldsymbol{u}(0))\Vert ,\end{eqnarray}$$

where the norms are Euclidean norms in $\mathbb{R}^{m}$ , and $\unicode[STIX]{x1D706}$ is defined as the LE corresponding to this CLV. CLVs with positive LEs are said to be unstable, CLVs with negative LEs are stable, and CLVs with zero LEs are neutral. The existence of CLVs is given by the Oseledets theorem (Oseledets Reference Oseledets1968), for which a proof in English can be found in Rezakhanlou (Reference Rezakhanlou2018, chap. 4), and more reference can be found at Oseledets (Reference Oseledets2008). In this paper, we sort LEs by decreasing order, so the $j$ th largest LE and its corresponding CLV will be called the $j$ th LE and $j$ th CLV, respectively.

Hyperbolicity is roughly defined as being that the tangent space at every state splits into a stable subspace, an unstable space and a neutral subspace. The one most widely used in ergodic theory is uniform hyperbolicity, which has two assumptions: (1) the three subspaces are uniformly bounded away from each other; and (2) the neutral subspace is one-dimensional. Our system cannot be uniform hyperbolic, since it violates the first assumption, that is, it has at least two neutral CLVs: the first one corresponds to the common time translation of continuous dynamical systems, and the second corresponds to spanwise translations due to the periodic boundary conditions.

Still, we want to check if our system violates the second assumption of uniform hyperbolicity. More specifically, we want to bound the norm of the operator which projects to one of the subspaces, i.e. $C_{\unicode[STIX]{x1D6FC}}$ given in Ni (Reference Ni2018b ), where we showed that $C_{\unicode[STIX]{x1D6FC}}=(1-\cos ^{2}\unicode[STIX]{x1D6FC})^{-1/2}$ , with $\unicode[STIX]{x1D6FC}$ being the smallest angle between subspaces. Of course, the ideal case would be that these subspaces are orthogonal to each other, leading to $\unicode[STIX]{x1D6FC}=90^{\circ }$ and $C_{\unicode[STIX]{x1D6FC}}=1$ takes its minimum. Both the bound of the norm of the shadowing direction (Pilyugin Reference Pilyugin1997) and the bound of the error in FD-NILSS (Ni & Wang Reference Ni and Wang2017) are proportional to $C_{\unicode[STIX]{x1D6FC}}$ , and we do not want it larger than 10, causing impact of an order of magnitude. The corresponding angle is $\unicode[STIX]{x1D6FC}_{0}=5.7^{\circ }$ , which serves as our threshold.

As we will see, the smallest angle between all CLVs can be affected by mesh or time length. However, we find a robust phenomenon related to hyperbolicity, that is, angles between two CLVs become much larger if their indices are further apart. We no longer have uniform hyperbolicity, since there is more than one neutral CLV and there may be tangencies between adjacent CLVs. On the other hand, this phenomenon is more robust, and reveals more structures of CLVs not indicated even by uniform hyperbolicity. In fact, we can even visualize this phenomenon by showing that active areas of CLVs are different.

More importantly, we will show that for our problem, where we do not have uniform hyperbolicity, tools given by ergodic theories are still valid, although they mostly use uniform hyperbolicity as logical assumptions. Specifically, we can still compute shadowing directions which correctly reflect sensitivities of statistics of the fluid system. Our results confirm the chaotic hypothesis (Gallavotti Reference Gallavotti2008), which says that many tools logically consequential to the uniform hyperbolicity assumption are in fact valid for a larger class of chaotic systems.

3.3 Numerical methods for computing Lyapunov exponents and covariant Lyapunov vectors

The algorithm we use to compute LEs is given by Benettin et al. (Reference Benettin, Galgani, Giorgilli and Strelcyn1980), and the algorithm for CLVs is by Ginelli et al. (Reference Ginelli, Poggi, Turchi, Chaté, Livi and Politi2007, Reference Ginelli, Chaté, Livi and Politi2013). Since these two algorithms share many of the same procedures, we describe them together. For $i=0,\ldots ,K-1$ , we define the $i$ th time segment as time span $[t_{i},t_{i+1}]$ , with $t_{i}=i\unicode[STIX]{x0394}T$ . In the algorithm presented below, for quantities defined on the entire segments such as $\boldsymbol{u}_{i}$ and $\unicode[STIX]{x1D652}_{i}$ , we use the same subscript as the segment they are defined on. For quantities defined only at interfaces between segments such as $\unicode[STIX]{x1D64C}_{i}$ and $\unicode[STIX]{x1D64D}_{i}$ , we use the same subscript as the time point they are defined at.

To begin, we should prescribe: (1) number of LEs/CLVs to compute, $M$ ; (2) length of each time segment, $\unicode[STIX]{x0394}T$ ; and (3) number of time segments, $K$ . Consequently, the time length of the entire trajectory, $T=K\unicode[STIX]{x0394}T$ , is determined. Then the algorithm is given by the following procedure.

1. 1. Generate initial conditions for primal solutions and homogeneous tangent solutions.

   1. 1.1. Compute the primal solution of (3.1) for sufficiently long time so that the trajectory lands onto the attractor, then set $t=0$ , and set initial condition of the primal system, $\boldsymbol{u}_{0}(0)$ .

   1. 1.2. Randomly generate an $m\times M$ orthogonal matrix $\unicode[STIX]{x1D64C}_{0}=[\boldsymbol{q}_{01},\ldots ,\boldsymbol{q}_{0M}]$ . This $\unicode[STIX]{x1D64C}_{0}$ will be used as initial condition for homogeneous tangent solutions.

2. 2. For $i=0$ to $K-1$ , on segment $i$ , where $t\in [t_{i},t_{i+1}]$ , do:

   1. 2.1. Compute the primal solution $\boldsymbol{u}_{i}(t)$ from $t_{i}$ to $t_{i+1}$ .

   2. 2.2. Compute homogeneous tangent solutions $\unicode[STIX]{x1D652}_{i}(t)=[\boldsymbol{w}_{i1}(t),\ldots ,\boldsymbol{w}_{iM}(t)]$ .

      1. 2.2.1. For each homogeneous tangent solution $\boldsymbol{w}_{ij}$ , $j=1,\ldots ,M$ , starting from initial condition $\boldsymbol{w}_{ij}(t_{i})=\boldsymbol{q}_{ij}$ , integrate (3.2) from $t_{i}$ to $t_{i+1}$ .

      2. 2.2.2. Perform QR factorization: $\unicode[STIX]{x1D652}_{i}(t_{i+1})=\unicode[STIX]{x1D64C}_{i+1}\unicode[STIX]{x1D64D}_{i+1}$ , where $\unicode[STIX]{x1D64C}_{i+1}=[\boldsymbol{q}_{i+1,1},\ldots ,\boldsymbol{q}_{i+1,M}]$ .

3. 3. The $j$ th largest LE, $\unicode[STIX]{x1D706}_{j}$ , is approximated by

   (3.5) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{j}=\frac{1}{K\unicode[STIX]{x0394}T}\mathop{\sum }_{i=1}^{K}\log |\unicode[STIX]{x1D60B}_{ij}|,\end{eqnarray}$$
   where $\unicode[STIX]{x1D60B}_{ij}$ is the $j$ th diagonal element in $\unicode[STIX]{x1D64D}_{i}$ . This formula for $\unicode[STIX]{x1D706}_{j}$ converges to the true value as $T$ becomes large.

4. 4. We define an $m\times M$ matrix $\unicode[STIX]{x1D651}(t)$ by

   (3.6) $$\begin{eqnarray}\unicode[STIX]{x1D651}(t)=\unicode[STIX]{x1D652}_{i}(t)\unicode[STIX]{x1D64D}_{i+1}^{-1}\cdots \unicode[STIX]{x1D64D}_{K}^{-1},\quad t\in [t_{i},t_{i+1}].\end{eqnarray}$$
   The $j$ th column of $\unicode[STIX]{x1D651}(t)$ converges to the direction of the $j$ th CLV when both $t$ and $(T-t)$ become large. Notice that, although $\unicode[STIX]{x1D651}(t)$ has different expressions on different segments, its columns are continuous across all segments.

We discuss the inner product we implicitly used in defining the orthogonality of the $\unicode[STIX]{x1D64C}$ matrices. The inner product between two flow fields at an instant is a summation of inner products between different components. For consistency of units, we normalize each component by free-stream properties. Moreover, we should prescribe a geometric metric for inner products between the same components. This metric should have bounded total volume, so that inner products are finite; it should also put more weight closer to the cylinder, so that FD-NILSS has better accuracy for the surface objectives we use later. A natural selection of such a metric is simply one where all mesh cells are equally weighted, and the entire space outside the mesh has weight zero. More specifically, with subscript $l$ denoting different components, and $k$ denoting different cells of the mesh, for $\boldsymbol{v},\boldsymbol{w}\in X=\mathbb{R}^{m}$ , we define inner products as

(3.7) $$\begin{eqnarray}\boldsymbol{v}^{\text{T}}\boldsymbol{w}:=\frac{1}{m^{\prime }}\mathop{\sum }_{k=1}^{m^{\prime }}\boldsymbol{v}_{k}^{\text{T}}\boldsymbol{w}_{k},\quad \text{where }\boldsymbol{v}_{k}^{\text{T}}\boldsymbol{w}_{k}:=\frac{1}{3}\left(\frac{v_{k1}w_{k1}}{\unicode[STIX]{x1D70C}_{0}^{2}}+\frac{\displaystyle \mathop{\sum }_{l=2}^{4}v_{kl}w_{kl}}{\unicode[STIX]{x1D70C}_{0}^{2}U_{0}^{2}}+\frac{v_{k5}w_{k5}}{\unicode[STIX]{x1D70C}_{0}^{2}E_{0}^{2}}\right).\end{eqnarray}$$

Here $m^{\prime }$ is the number of cells. We order the five components of states by $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D70C}\boldsymbol{U}$ , $\unicode[STIX]{x1D70C}E$ , and normalize all components of $\unicode[STIX]{x1D70C}\boldsymbol{U}$ by $\unicode[STIX]{x1D70C}U_{0}$ . Under our definition, if every point in a flow field has the same state as the free stream, then this flow field has unit norm.

The convergence of CLV in the above algorithm is in terms of the distance between one-dimensional subspaces. Since a CLV is a homogeneous tangent solution, we can multiply it by any factor, and still get a homogeneous tangent solution whose norm behaves like an exponential function with the same LE, which, by definition, is still the same CLV. In other words, it is only the directions of CLVs that are meaningful, but not the magnitudes. Hence we can normalize CLVs by any factor we like. In this paper, we normalize CLVs such that their maximal values are 1.0.

In this paper, we use finite differences to approximate tangent solutions. To approximate a homogeneous solution $w$ with initial condition $\boldsymbol{w}^{0}$ , we use the definition $\boldsymbol{w}=\text{d}\boldsymbol{u}/\text{d}\unicode[STIX]{x1D719}\approx \unicode[STIX]{x0394}\boldsymbol{u}/\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ , where $\unicode[STIX]{x1D719}$ controls the perturbation in the initial condition in the direction of $\boldsymbol{w}^{0}$ . More specifically, we compute a perturbed primal solution $\boldsymbol{u}^{\boldsymbol{w}}$ with a perturbed initial condition $\boldsymbol{u}^{0}+\boldsymbol{w}^{0}\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ , and the approximation for $\boldsymbol{w}$ is $\boldsymbol{w}\approx (\boldsymbol{u}^{\boldsymbol{w}}-\boldsymbol{u})/\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ , where $\boldsymbol{u}$ is the base trajectory.

3.4 Results of Lyapunov exponents

The time-step size used for computing LEs and CLVs is the same as in the numerical simulation, that is, $\unicode[STIX]{x0394}t=10^{-8}=1.32\times 10^{-3}t_{0}$ . The number of homogeneous tangent solutions we compute is $M=40$ . Each segment has 200 time steps, which gives $\unicode[STIX]{x0394}T=2\times 10^{-6}=0.264t_{0}$ . The total number of time segments is $K=800$ , hence the time length of the entire trajectory is $T=1.6\times 10^{-3}=211t_{0}$ .

With the above numerical settings, the convergence history of the first 40 LEs is shown on the left of figure 5. As explained in Ni & Wang (Reference Ni and Wang2017), the confidence interval of an LE is estimated by the smallest interval which bounds the history of the LE and whose size shrinks as $T^{-0.5}$ . The first 40 LEs and their confidence intervals are shown on the right of figure 5.

Figure 5. Lyapunov exponents (LE) normalized by $t_{0}^{-1}$ . Left: convergence history of the first 40 different LEs. Right: confidence intervals of LEs. (a) Results on the coarser mesh. The largest LE is $0.20t_{0}^{-1}$ , meaning in one flow-through time $t_{0}$ , the norm of the first CLV becomes $\text{e}^{0.20}=1.22$ times larger. (b) Results on the finer mesh. The largest LE is $0.21t_{0}^{-1}$ .

We can see that LE results depend on the meshes, as similarly discussed by Fernandez & Wang (Reference Fernandez and Wang2017b ). In figure 5, although the general shape of the LE spectrum looks similar, a closer look will show that the finer mesh has larger but fewer positive LEs. As a result, the selection of neutral CLVs is ambiguous and changes due to meshes.

However, we will show that there are several important physical phenomena which, although they are logically consequential to the Lyapunov analysis results, exist robustly regardless of the detailed value of the LE spectrum. This agrees with the chaotic hypothesis (Gallavotti Reference Gallavotti2008), which says that many tools logically consequential to the uniform hyperbolicity assumption are in fact valid for a larger class of dynamical systems. More specifically, those phenomena are as follows:

1. (1) the dimension of the attractor is much lower than the dimension of the system;

2. (2) angles between far-apart CLVs are large, and, in particular, CLVs are active at different locations;

3. (3) shadowing directions exist, and reflect the sensitivities of statistics to parameter perturbations.

The focus of this paper will be the last two phenomena. Moreover, for engineering interests, we will show that FD-NILSS can compute shadowing directions and, further, the sensitivities with reasonable cost.

We first show that the dimension of the attractor is much lower than the dimension of the system. Intuitively, to remain its invariant hyper-volume, the attractor manifold should contain enough contracting (stable) CLVs to balance out the expansion caused by unstable CLVs. More specifically, Frederickson et al. (Reference Frederickson, Kaplan, Yorke and Yorke1983) estimated the attractor dimension by the so-called Lyapunov dimension:

(3.8) $$\begin{eqnarray}D_{\unicode[STIX]{x1D706}}=N+\frac{1}{|\unicode[STIX]{x1D706}_{N+1}|}\mathop{\sum }_{j=1}^{j=N}\unicode[STIX]{x1D706}_{j}\in [N,N+1],\end{eqnarray}$$

where $N$ is such that $\unicode[STIX]{x1D706}_{1}+\unicode[STIX]{x1D706}_{2}+\cdots +\unicode[STIX]{x1D706}_{N}>0$ and $\unicode[STIX]{x1D706}_{1}+\unicode[STIX]{x1D706}_{2}+\cdots +\unicode[STIX]{x1D706}_{N+1}<0$ . We want to give a bound to $N$ on the finer mesh. To do this, we remind readers that the algorithm we are using computes LEs by descending order, so we have $\unicode[STIX]{x1D706}_{j}\leqslant \unicode[STIX]{x1D706}_{40}$ for all $j\geqslant 40$ . Using the first condition of $N$ , we have

(3.9) $$\begin{eqnarray}\mathop{\sum }_{j=1}^{40}\unicode[STIX]{x1D706}_{j}+(N-40)\unicode[STIX]{x1D706}_{40}=\mathop{\sum }_{j=1}^{40}\unicode[STIX]{x1D706}_{j}+\mathop{\sum }_{j=41}^{N}\unicode[STIX]{x1D706}_{40}\geqslant \mathop{\sum }_{j=1}^{40}\unicode[STIX]{x1D706}_{j}+\mathop{\sum }_{j=41}^{N}\unicode[STIX]{x1D706}_{j}>0,\end{eqnarray}$$

which further leads to

(3.10) $$\begin{eqnarray}N\leqslant 40+\frac{\displaystyle \mathop{\sum }_{j=1}^{40}\unicode[STIX]{x1D706}_{j}}{|\unicode[STIX]{x1D706}_{40}|}=40+\frac{1.031}{0.027}=78.\end{eqnarray}$$

On the other hand, since $\unicode[STIX]{x1D706}_{1}+\cdots +\unicode[STIX]{x1D706}_{40}>0$ , by the second property of $N$ , we know $N>40$ . Hence the Lyapunov dimension $41\leqslant D_{\unicode[STIX]{x1D706}}\leqslant 79$ . In other words, the chaotic dynamic of our flow problem on the finer mesh can be attributed to the interaction of fewer than 79 degrees of freedom. Using the same method, we can estimate that the attractor dimension on the coarser mesh is smaller than 109.

3.5 Results of covariant Lyapunov vectors and hyperbolicity

In figure 6(a,c), we plot the histogram of angles between all pairs of the first 40 CLVs, which are perturbations whose norm grows exponentially, as defined in (3.4). Notice the range of angles is $[0^{\circ },90^{\circ }]$ , since the angle between the directions of two CLVs is the angle between subspaces of dimension one. For the finer mesh, the smallest angle is $16.4^{\circ }$ . However, for the coarser mesh, the smallest angle is only $5.8^{\circ }$ , almost equal the threshold value of $5.7^{\circ }$ . This indicates that the smallest angle and hence the second assumption in uniform hyperbolicity is not a robust criterion to check for fluid systems, since it is plausible that the smallest angle falls below the threshold for another mesh or a longer time span.

Figure 6. Histogram of angles between CLVs on segments 350–449, with the entire trajectory of length $T=211t_{0}$ partitioned into $K=800$ segments. The total area of the histogram is normalized to one. (a) Coarser mesh, all CLVs, smallest angle is $5.8^{\circ }$ . (b) Coarser mesh, CLV indices more than five apart, smallest angle is $29.2^{\circ }$ . (c) Finer mesh, all CLVs, smallest angle is $16.4^{\circ }$ . (d) Finer mesh, CLV indices more than five apart, smallest angle is $40.8^{\circ }$ .

The above dependences on meshes call for a closer examination to find some robust physical phenomena not depending on meshes. A deeper look into CLV results finds that tangencies happen only among adjacent CLVs: this result is similar to that obtained by Xu & Paul (Reference Xu and Paul2016). Hence we plot the angles between CLVs whose indices are more than five apart. For example, for the 20th CLV, we consider its angles with the 1st to 15th CLVs and the 25th to 40th CLVs. This result is plotted in figure 6(b,d).

As we can see, when the exponents of CLVs are further apart, their angles become larger: this phenomenon is related to hyperbolicity but not formally named yet. So, on one hand, what we observe violates the second assumption of uniform hyperbolicity since occasional tangencies may fail the threshold. But, on the other hand, our observation is even stronger than the uniform hyperbolicity, because it specifies a trend not suggested by the definition of uniform hyperbolicity.

To further illustrate, we plot minimal and averaged angles between all pairs of CLVs in figure 7. Not surprisingly, most tangencies happen between adjacent CLVs, and the angles between CLVs get larger for CLVs whose indices are farther away. In particular, on the coarser mesh, the first seven CLVs are closely coupled, whereas the first several CLVs on the finer mesh are less tangent to each other. In general, tangencies among CLVs are rare events for both meshes, but even less frequent for the finer mesh. Once again, the rareness of the tangencies indicates that smallest angle and hence uniform hyperbolicity is not a robust criterion; on the other hand, the fact that angles between far-apart CLVs are large is robust to meshes, and this may lead to some definition of a new form of hyperbolicity.

Figure 7. Minimal and averaged angles between all pairs of CLVs. The minimum and average are taken on segments 350–449, with the entire trajectory of length $T=211t_{0}$ partitioned into $K=800$ segments. (a) Coarser mesh, minimal angles between all pairs of CLVs. (b) Coarser mesh, averaged angles between all pairs of CLVs. (c) Finer mesh, minimal angles between all pairs of CLVs. (d) Finer mesh, averaged angles between all pairs of CLVs.

In fact, for two CLVs whose exponents are apart, their angles are so large that we can even observe the difference in active areas. To illustrate, for the coarser mesh, we plot the 1st and 40th CLVs in figure 8. For the finer mesh, we plot the 1st, 5th, 17th and 40th CLVs in figure 9. Movies of the 1st, 5th, 17th and 28th CLVs on the finer mesh are movie 1 to movie 4 of the supplementary material, available at https://doi.org/10.1017/jfm.2018.986. We remind readers that only the direction of CLVs, that is, only the relative flow pattern of CLVs, is meaningful. Since our primal solution lives in the function space $X\subset \mathbb{R}^{m}$ whose tangent space is also $\mathbb{R}^{m}$ , for any time $t$ , a CLV $\unicode[STIX]{x1D73B}(t)$ lives in a function space locally the same as the primal solution $\boldsymbol{u}(t)$ . In our 3D flow problem, this means our CLVs also have $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D70C}\boldsymbol{U}$ , $\unicode[STIX]{x1D70C}E$ components, and we plot $|\unicode[STIX]{x1D70C}\boldsymbol{U}|$ .

Figure 8. CLVs computed on the coarser mesh on a trajectory of time length $T=211t_{0}$ : (a) 1st CLV; (b) 40th CLV; and (c) primal flow field. Left: snapshots at $T/2$ . Right: averaged over time span [ $95t_{0},116t_{0}$ ]. Plotted by the magnitude of $\unicode[STIX]{x1D70C}\boldsymbol{U}$ component, and normalized such that the largest value is 1.00. The primal flow at the bottom is normalized by $\unicode[STIX]{x1D70C}_{0}U_{0}$ .

Figure 9. CLV computed on the finer mesh, using same settings as the coarser mesh: (a) 1st CLV; (b) 5th CLV; (c) 17th CLV; (d) 40th CLV; and (e) primal flow field.

As we can see in both figures 8 and 9, the general trend is that the active areas of CLVs move downstream as exponents become smaller. Intuitively, in boundary layers and near wakes, chaos is generated; perturbations in these places tend to grow, resulting in unstable CLVs. On the other hand, in far wakes, chaos is dissipated; perturbations in this region tend to decay, resulting in stable CLVs. Further, we can rationally conjecture that the location of neutral CLVs is in between stable and unstable CLVs, that is, their active area is the entire wake. This conjecture is also backed by observing the plot of a neutral CLV in figure 9(c).

We conjecture that perturbations in the stagnant area in front of the cylinder will be dissipated, resulting in stable CLVs. This is because the fluid motion in front of the cylinder is almost determined, that is, fluid here should decelerate into stagnancy. If we make a perturbation to this area, the fluid motion in the long term will not be changed. Hence, such a perturbation could not last long, and the corresponding LE should be negative. Our conjecture is also partly evidenced by the fact that the first 40 CLVs are not active in this area.

We further conjecture that perturbations in the free stream are also dissipated, resulting in stable CLVs. We provide two arguments supporting this conjecture. First, assume we make a perturbation in the free stream. Then, after changing reference frame, this is equivalent to making a small perturbation in a stagnant fluid domain; the Euclidean norm of such a perturbation decays due to the dissipation in the Navier–Stokes equations. Second, since our metric, induced by (3.7), is compactly supported, a perturbation in the free stream will be transported outside the mesh where the volume weight is zero. Hence under our norm, perturbation in the free stream decays also because of the decay in the volume weight. In fact, because our norm of such perturbation decays to zero, the corresponding LE should be very negative.

To summarize, we make the following conjecture about active areas of CLVs in subsonic open flows. In subsonic open flows, unstable CLVs are active in the instability-generating area such as boundary layers and near wakes; the neutral CLV occupies similar locations as the wake; stable CLVs with moderately negative LEs are active in the dissipative area such as far wakes; stable CLVs active in the stagnant area ahead of the cylinder and in the free stream have very negative LEs. Our conjecture also indicates that subsonic open flows have large angles between far-apart CLVs.

Our conjecture may provide insights into the question about how many CLVs are unstable for fluid systems. Owing to our arguments, for open flows, or more generally fluid systems with non-turbulent regions, there should be a large part of CLVs that are stable; hence unstable CLVs should be a small fraction of all CLVs. At least for these fluid systems, more analysis should be focused on unstable CLVs since they are the major contributor for turbulence. Additionally, there should be more numerical methods developed which explicitly exploit the unstable structures, such as NILSS and NILSAS.

We find more structures in hyperbolicity, but our system certainly does not satisfy the uniform hyperbolicity under which shadowing methods were mathematically developed. Yet there is still hope, as the chaotic hypothesis suggested that, although the logical bridge of uniform hyperbolicity breaks down, we may find tools across the river still useful. With this hope, we go on to examine shadowing results on our fluid system. Should we find shadowing methods still valid, we will have more confidence in its generality.

4.1 Definition of shadowing directions

The shadowing solution $\boldsymbol{v}^{\infty }$ is an inhomogeneous tangent solution whose orthogonal projection perpendicular to the trajectory, $\boldsymbol{v}^{\infty \bot }$ , is uniformly bounded on an infinitely long trajectory. Here $(\cdot )^{\bot }$ is defined as

(4.1) $$\begin{eqnarray}\boldsymbol{p}^{\bot }=\boldsymbol{p}-\frac{\boldsymbol{f}^{\text{T}}\boldsymbol{p}}{\boldsymbol{f}^{\text{T}}\boldsymbol{f}}\,\boldsymbol{f},\end{eqnarray}$$

where $\boldsymbol{p}\in \mathbb{R}^{m}$ is an arbitrary vector, $\boldsymbol{f}$ is the trajectory direction defined in (3.1), and the inner product was defined in (3.7). We remark here that, for a continuous dynamical system, the definition of shadowing direction involves a reparametrization of time (Pilyugin Reference Pilyugin1999, Definition 1.26), or equivalently a ‘time dilation’ (Wang et al. Reference Wang, Hu and Blonigan2014), or projecting the tangent solutions onto the subspace perpendicular to $\boldsymbol{f}$ (Ni & Wang Reference Ni and Wang2017). The $\boldsymbol{v}^{\infty \bot }$ is the projection of $\boldsymbol{v}^{\infty }$ given by the third approach, and $\boldsymbol{v}^{\infty \bot }$ is bounded, but $\boldsymbol{v}^{\infty }$ is not. For convenience, we also call $\boldsymbol{v}^{\bot }$ the shadowing direction.

The shadowing direction has the following physical meaning. The definition of inhomogeneous tangent solution $\boldsymbol{v}=\text{d}\boldsymbol{u}/\text{d}s$ tells us that, with perturbed parameter $s+\unicode[STIX]{x1D6FF}s$ , there is a new trajectory $\boldsymbol{u}+\unicode[STIX]{x1D6FF}\boldsymbol{u}$ such that $\unicode[STIX]{x1D6FF}\boldsymbol{u}\approx \boldsymbol{v}^{\infty }\unicode[STIX]{x1D6FF}s$ . Now the perpendicular distance between the new and the old trajectory is $\unicode[STIX]{x1D6FF}\boldsymbol{u}^{\bot }\approx \boldsymbol{v}^{\infty \bot }\unicode[STIX]{x1D6FF}s$ ; together with the boundedness of $\boldsymbol{v}^{\infty \bot }$ , we see that the new trajectory remains perpendicularly close to the old trajectory for a long time. This new trajectory will be referred to as the shadowing trajectory, and the old one as the base trajectory. This intuition of the shadowing direction is shown in figure 10.

The initial condition of the shadowing direction is not known a priori, or equivalently, we do not specify how system parameters would affect initial conditions. In fact, as we shall see, although the shadowing direction solves that same ODE as in (3.3), its determining condition is no longer the initial condition; rather, to determine a shadowing direction, we minimize the $L^{2}$ norm of the orthogonal projection of an inhomogeneous tangent solution.

Figure 10. Shadowing directions. The base trajectory has parameter $s$ , the shadowing trajectory has parameter $s+\unicode[STIX]{x1D6FF}s$ . The first-order approximation of the perpendicular distance is $\unicode[STIX]{x1D6FF}\boldsymbol{u}^{\bot }\approx \boldsymbol{v}^{\infty \bot }\unicode[STIX]{x1D6FF}s$ , where $\boldsymbol{v}^{\infty \bot }$ is the shadowing direction.

4.2 Sensitivity analysis of long-time averages via shadowing directions

For chaotic dynamical systems such as the 3D flow problem in this paper, the output of the system, such as the drag or lift, is typically aperiodic. As often happens in engineering, the objective is the long-time average $\langle J\rangle _{\infty }$ of an instantaneous quantity $J(\boldsymbol{u},s):\mathbb{R}^{m}\times \mathbb{R}\rightarrow \mathbb{R}$ . More specifically, we define

(4.2) $$\begin{eqnarray}\langle J\rangle _{\infty }:=\mathop{\lim }_{T\rightarrow \infty }\langle J\rangle _{T},\quad \text{where }\langle J\rangle _{T}:=\frac{1}{T}\int _{0}^{T}J(\boldsymbol{u},s)\,\text{d}t.\end{eqnarray}$$

Typically $\langle J\rangle _{\infty }$ is approximated by a $\langle J\rangle _{T}$ with a large $T$ . The sensitivity of the objective with respect to system parameters, $\text{d}\langle J\rangle _{\infty }/\text{d}s$ , is of engineering interest. In this subsection we show how to compute the sensitivity using the shadowing direction.

Assume we have found, on a finitely long trajectory, an approximation $\boldsymbol{v}\approx \boldsymbol{v}^{\infty }$ , where $\boldsymbol{v}$ is also an inhomogeneous tangent solution. We first define the time dilation term (Wang et al. Reference Wang, Hu and Blonigan2014) $\unicode[STIX]{x1D702}$ :

(4.3) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{v}^{\bot }}{\text{d}t}=\unicode[STIX]{x2202}_{\boldsymbol{u}}\,f\boldsymbol{v}^{\bot }+\unicode[STIX]{x2202}_{s}\,\boldsymbol{f}+\unicode[STIX]{x1D702}\boldsymbol{f}.\end{eqnarray}$$

Intuitively, $\unicode[STIX]{x1D702}$ describes the relative time spent on the shadowing trajectory in comparison with the base trajectory. In figure 10, if the shadowing trajectory takes less time to travel from point C to D than the base trajectory from A to B, then (4.3) gives $\unicode[STIX]{x1D702}<0$ . On the other hand, if the shadowing trajectory moves slower, then $\unicode[STIX]{x1D702}>0$ .

Once we get $\boldsymbol{v}$ and $\unicode[STIX]{x1D702}$ , the sensitivity can be computed by

(4.4) $$\begin{eqnarray}\frac{\text{d}\langle J\rangle _{\infty }}{\text{d}s}\approx \frac{1}{T}\int _{0}^{T}[\unicode[STIX]{x2202}_{\boldsymbol{ u}}J\boldsymbol{v}^{\bot }+\unicode[STIX]{x2202}_{s}J+\unicode[STIX]{x1D702}(J-\langle J\rangle _{T})]\,\text{d}t,\end{eqnarray}$$

where $\langle J\rangle _{T}$ is defined in (4.2). Intuitively, the first term inside the integration in (4.4) describes the contribution by the perpendicular distance between the shadowing trajectory and the base trajectory. The second term is due to the fact that the function $J$ may explicitly depend on $s$ . The last term involves $\unicode[STIX]{x1D702}$ , and it accounts for the fact that the average is taken with respect to time: if in some small region the shadowing trajectory moves faster than the base trajectory, then the time spent is shorter, hence the contribution from this region in the time average should also be smaller.

Another formula for the sensitivity is easier for computer programming:

(4.5) $$\begin{eqnarray}\frac{\text{d}\langle J\rangle _{\infty }}{\text{d}s}\approx \frac{1}{T}\left[\int _{0}^{T}(\unicode[STIX]{x2202}_{\boldsymbol{ u}}J\boldsymbol{v}+\unicode[STIX]{x2202}_{s}J)\,\text{d}t+\unicode[STIX]{x1D709}|_{0}^{T}\langle J\rangle _{T}-(\unicode[STIX]{x1D709}J)|_{0}^{T}\right],\end{eqnarray}$$

where the time difference term, $\unicode[STIX]{x1D709}$ , is a time-dependent scalar function such that

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D709}\boldsymbol{f}=\boldsymbol{v}-\boldsymbol{v}^{\bot }.\end{eqnarray}$$

Thus $\unicode[STIX]{x1D709}$ is easier to compute than $\unicode[STIX]{x1D702}$ , since its definition does not involve a time derivative. Note that in (4.5) we use $\boldsymbol{v}$ instead of its projection $\boldsymbol{v}^{\bot }$ . The derivation of (4.4), (4.5) and their relation can be found in the appendix of Ni & Wang (Reference Ni and Wang2017).

4.3 The NILSS algorithm for computing shadowing directions

We describe the non-intrusive least-squares shadowing (NILSS) algorithm in this subsection. As in § 3.3, for $i=0,\ldots ,K-1$ , we define the $i$ th time segment as $[t_{i},t_{i+1}]$ , with $t_{i}=i\unicode[STIX]{x0394}T$ . In the algorithm presented below, for quantities defined on the entire segments, such as $\boldsymbol{u}_{i},~\unicode[STIX]{x1D652}_{i},~\boldsymbol{v}_{i}^{\ast }$ , $\unicode[STIX]{x1D63E}_{i}$ and $\boldsymbol{d}_{i}$ , we use the same subscripts as the segments they are defined on. For quantities defined only at interfaces between segments, such as $\unicode[STIX]{x1D64C}_{i},~\unicode[STIX]{x1D64D}_{i}$ and $\boldsymbol{b}_{i}$ , we use the same subscripts as the time points they are defined at. Again, we use finite difference results to approximate tangent solutions: such a variant is called the finite difference non-intrusive least-squares shadowing (FD-NILSS) algorithm, whose details are given in Ni et al. (Reference Ni, Wang, Fernandez and Talnikar2018).

To start with, we should prescribe: (1) number of homogeneous tangent solutions, $M$ , which must be larger than the number of unstable CLVs, $m_{us}$ ; (2) length of each time segment, $\unicode[STIX]{x0394}T$ ; and (3) number of time segments, $K$ . Consequently, the time length of the entire trajectory, $T=K\unicode[STIX]{x0394}T$ , is also determined. Then, the NILSS algorithm is given by the following procedure.

1. 1. Generate initial conditions for primal solutions, homogeneous tangent solutions and inhomogeneous tangent solutions.

   1. 1.1. Compute the primal solution of (3.1) for sufficiently long time so that the trajectory lands on the attractor, then set $t=0$ , and set the initial condition of the primal system, $\boldsymbol{u}_{0}(0)$ .

   1. 1.2. Randomly generate an $m\times M$ orthogonal matrix $\unicode[STIX]{x1D64C}_{0}=[\boldsymbol{q}_{01},\ldots ,\boldsymbol{q}_{0M}]$ whose column vectors are orthogonal to $\boldsymbol{f}(t=0)$ . This $\unicode[STIX]{x1D64C}_{0}$ will be used as initial conditions for homogeneous tangent solutions.

   1. 1.3. Set the initial condition of the particular inhomogeneous tangent solution $\boldsymbol{v}_{0}^{\ast }(0)=0$ .

2. 2. For $i=0$ to $K-1$ , on segment $i$ , where $t\in [t_{i},t_{i+1}]$ , do:

   1. 2.1. Compute the primal solution $\boldsymbol{u}_{i}(t)$ from $t_{i}$ to $t_{i+1}$ .

   2. 2.2. Compute homogeneous tangent solutions $\unicode[STIX]{x1D652}_{i}(t)=[\boldsymbol{w}_{i1}(t),\ldots ,\boldsymbol{w}_{iM}(t)]$ .

      1. 2.2.1. For each homogeneous tangent solution $\boldsymbol{w}_{ij}$ , $j=1,\ldots ,M$ , starting from initial condition $\boldsymbol{w}_{ij}(t_{i})=\boldsymbol{q}_{ij}$ , integrate (3.2) from $t_{i}$ to $t_{i+1}$ .

      2. 2.2.2. Compute orthogonal projection $\unicode[STIX]{x1D652}_{i}^{\bot }(t)=[\boldsymbol{w}_{i1}^{\bot }(t),\ldots ,\boldsymbol{w}_{iM}^{\bot }(t)]$ via (4.1).

      3. 2.2.3. Compute and store the covariant matrix on segment $i$ :

         (4.7) $$\begin{eqnarray}\unicode[STIX]{x1D63E}_{i}=\int _{t_{i}}^{t_{i+1}}(\unicode[STIX]{x1D652}_{i}^{\bot })^{\text{T}}\unicode[STIX]{x1D652}_{i}^{\bot }\,\text{d}t.\end{eqnarray}$$

      4. 2.2.4. Perform QR factorization: $\unicode[STIX]{x1D652}_{i}^{\bot }(t_{i+1})=\unicode[STIX]{x1D64C}_{i+1}\unicode[STIX]{x1D64D}_{i+1}$ , where $\unicode[STIX]{x1D64C}_{i+1}=[\boldsymbol{q}_{i+1,1},\ldots ,\boldsymbol{q}_{i+1,M}]$ .

   3. 2.3. Compute the particular inhomogeneous tangent solution $\boldsymbol{v}_{i}^{\ast }(t)$ .

      1. 2.3.1. Starting from initial condition $\boldsymbol{v}_{i}^{\ast }(t_{i})$ , integrate the inhomogeneous equation (3.3) from $t_{i}$ to $t_{i+1}$ .

      2. 2.3.2. Compute the orthogonal projection $\boldsymbol{v}_{i}^{\ast \bot }(t)$ via (4.1).

      3. 2.3.3. Compute and store

         (4.8) $$\begin{eqnarray}\boldsymbol{d}_{i}=\int _{t_{i}}^{t_{i+1}}(\unicode[STIX]{x1D652}_{i}^{\bot })^{\text{T}}\boldsymbol{v}_{i}^{\ast \bot }\,\text{d}t.\end{eqnarray}$$

      4. 2.3.4. Orthogonalize $\boldsymbol{v}_{i}^{\ast \bot }(t_{i+1})$ with respect to $\unicode[STIX]{x1D652}_{i+1}^{\bot }(t_{i+1})=\unicode[STIX]{x1D64C}_{i+1}$ to obtain the initial condition of the next time segment:

         (4.9) $$\begin{eqnarray}\boldsymbol{v}_{i+1}^{\ast }(t_{i+1})=\boldsymbol{v}_{i}^{\ast \bot }(t_{i+1})-\unicode[STIX]{x1D64C}_{i+1}\boldsymbol{b}_{i+1},\end{eqnarray}$$
         compute and store:
         (4.10) $$\begin{eqnarray}\boldsymbol{b}_{i+1}=\unicode[STIX]{x1D64C}_{i+1}^{\text{T}}\boldsymbol{v}_{i}^{\ast \bot }(t_{i+1}).\end{eqnarray}$$

3. 3. Solve the NILSS problem:

   (4.11) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\mathop{\min }_{\{a_{i}\}}\displaystyle \mathop{\sum }_{i=0}^{K-1}\frac{1}{2}\boldsymbol{a}_{i}^{\text{T}}\unicode[STIX]{x1D63E}_{i}\boldsymbol{a}_{i}+\boldsymbol{d}_{i}^{\text{T}}\boldsymbol{a}_{i}\\ \text{s.t.}\quad \boldsymbol{a}_{i}=\unicode[STIX]{x1D64D}_{i}\boldsymbol{a}_{i-1}+\boldsymbol{b}_{i},\quad i=1,\ldots ,K-1.\end{array}\right\}\end{eqnarray}$$
   This is a least-squares problem with arguments $\{\boldsymbol{a}_{i}\}\text{}_{i=0}^{i=K-1}$ , where $\boldsymbol{a}_{i}\in \mathbb{R}^{M}$ for each  $i$ .

4. 4. Compute $v_{i}$ within each time segment $t\in [t_{i},t_{i+1}]$ :

   (4.12) $$\begin{eqnarray}\boldsymbol{v}_{i}(t)=\boldsymbol{v}_{i}^{\ast }(t)+\unicode[STIX]{x1D652}_{i}(t)\boldsymbol{a}_{i}.\end{eqnarray}$$
   The orthogonal projections of $\boldsymbol{v}_{i}$ on each segment, $\{\boldsymbol{v}_{i}^{\bot }\}\text{}_{i=0}^{K-1}$ , c onverge to the orthogonal projection of the shadowing direction, $\boldsymbol{v}^{\infty \bot }$ , as both $t$ and $T-t$ become large. In other words, we have smaller error in the middle section of a long time span.

5. 5. Compute $\unicode[STIX]{x1D709}_{i}$ at the end of each segment:

   (4.13) $$\begin{eqnarray}\unicode[STIX]{x1D709}_{i}(t_{i+1})=\frac{(\boldsymbol{v}_{i}(t_{i+1}))^{\text{T}}\boldsymbol{f}(\boldsymbol{u}(t_{i+1}))}{\boldsymbol{f}(\boldsymbol{u}(t_{i+1}))^{\text{T}}\boldsymbol{f}(\boldsymbol{u}(t_{i+1}))}.\end{eqnarray}$$

6. 6. The derivative can be computed by

   (4.14) $$\begin{eqnarray}\frac{\text{d}\langle J\rangle _{\infty }}{\text{d}s}\approx \frac{1}{T}\mathop{\sum }_{i=0}^{K-1}\left[\int _{t_{i}}^{t_{i+1}}(\unicode[STIX]{x2202}_{\boldsymbol{u}}J\boldsymbol{v}_{i}+\unicode[STIX]{x2202}_{s}J)\,\text{d}t+\unicode[STIX]{x1D709}_{i}(t_{i+1})(\langle J\rangle _{T}-J(t_{i+1}))\right].\end{eqnarray}$$

Readers should note that, although at the beginning of each segment, $\boldsymbol{v}_{i}$ is perpendicular to $\boldsymbol{f}(t_{i})$ , typically $\boldsymbol{v}_{i}$ is not perpendicular to $\boldsymbol{f}(t_{i})$ for the rest of the segment. In NILSS, we typically use perpendicularly projected vectors, such as $\boldsymbol{v}_{i}^{\ast \bot }$ , $\unicode[STIX]{x1D652}_{i}^{\bot }$ , $\boldsymbol{v}_{i}^{\bot }$ , for constructing the minimization problem in (4.11) to compute the shadowing directions. Unprojected vectors, such as $\boldsymbol{v}_{i}$ , are typically used for computing sensitivities by (4.14). As shown in the appendix of Ni & Wang (Reference Ni and Wang2017), we can construct a continuous $\boldsymbol{v}$ , and its orthogonal projection $\boldsymbol{v}^{\bot }$ confined on each segment yields $\boldsymbol{v}_{i}^{\bot }$ . However, $\boldsymbol{v}_{i}$ are not confinements of $\boldsymbol{v}$ , and $\boldsymbol{v}_{i}$ are not necessarily continuous across segments. In fact, in (4.11), the minimization is minimizing $\Vert \boldsymbol{v}^{\bot }\Vert _{L^{2}}$ , while the constraint prescribes that $\boldsymbol{v}^{\bot }$ is continuous across all time segments.

Intuitively, at $t=0$ , the difference between $\boldsymbol{v}^{\infty \bot }$ and $\boldsymbol{v}^{\bot }$ contains unstable CLVs; since $\boldsymbol{v}^{\infty \bot }$ is bounded, as time evolves, the exponential growth of unstable CLVs increases both $\Vert \boldsymbol{v}^{\infty \bot }-\boldsymbol{v}^{\bot }\Vert _{L^{2}}$ and $\Vert \boldsymbol{v}^{\ast \bot }\Vert _{L^{2}}$ . Also, as time evolves, the span of $M$ homogeneous tangent solutions approximates the span of the first $M$ CLVs, hence $\boldsymbol{v}^{\bot }=\boldsymbol{v}^{\ast \bot }+\unicode[STIX]{x1D652}^{\bot }\boldsymbol{a}$ allows us to subtract unstable CLVs from $\boldsymbol{v}^{\bot }$ . Now minimizing $\Vert \boldsymbol{v}^{\bot }\Vert _{L^{2}}$ removes unstable CLVs from both $\boldsymbol{v}^{\bot }$ and $\boldsymbol{v}^{\infty \bot }-\boldsymbol{v}^{\bot }$ , thus making $\boldsymbol{v}_{i}^{\bot }\approx \boldsymbol{v}^{\infty \bot }$ . Since the number of unstable CLVs is typically much smaller than the dimension of the dynamical system, NILSS is computationally efficient. A more detailed explanation of NILSS is contained in Ni & Wang (Reference Ni and Wang2017).

4.4 Results of shadowing directions

This subsection discusses the shadowing directions computed by the FD-NILSS algorithm. We remind readers that shadowing solutions depend only on the choice of system parameters but not on objectives. In particular, this subsection shows two shadowing solutions with respect to two system parameters: (1)  $x$ -directional free-stream velocity $U$ , normalized by $U_{0}$ ; and (2) rotation speed of the cylinder $\unicode[STIX]{x1D714}$ , measured in cycles per time unit, normalized by  $\unicode[STIX]{x1D714}_{0}$ .

In FD-NILSS we set time-step size $1\times 10^{-8}$ , and 200 steps in each time segment; hence the segment length is $\unicode[STIX]{x0394}T=2\times 10^{-6}=0.259t_{0}$ . These are the same as what we used for computing LEs and CLVs. Different from those used for computing LEs and CLVs, we set the number of homogeneous tangents $M=30$ and the number of segments $K=600$ ; hence the time length of the entire trajectory is $T=1.2\times 10^{-3}=158.4t_{0}$ .

We first plot the norm of shadowing directions on the finer mesh in figure 11, where the norm is induced by the inner product defined in (3.7). Since $\boldsymbol{v}^{\bot }=\text{d}\boldsymbol{u}^{\bot }/\text{d}s$ , with $\text{d}\boldsymbol{u}^{\bot }$ already normalized by free-stream constants in the definition of inner product, we further need to normalize shadowing directions by $s_{0}^{-1}$ . Here $s_{0}$ is the unit parameter, that is, $U_{0}$ or $\unicode[STIX]{x1D714}_{0}$ .

Figure 11. Norms of shadowing directions on the finer mesh, measured from $66t_{0}$ to $93t_{0}$ , computed on a trajectory of length $T=158t_{0}$ : (a) with $U$ as parameter; and (b) with $\unicode[STIX]{x1D714}$ as parameter.

By inspecting figure 11, we can verify the uniform boundedness of shadowing directions, as both $\Vert \boldsymbol{v}^{\bot }\Vert$ do not change much from $66t_{0}$ to $93t_{0}$ . For comparison we may think of the first CLV, whose norm would grow $\text{e}^{0.21\times (93-66)}=290$ times larger within the same time span.

Orthogonal projections of shadowing directions, $\boldsymbol{v}^{\bot }$ , are plotted in figure 12. Movies of $\boldsymbol{v}^{\bot }$ on the finer mesh are movie 5 and movie 6 of the supplementary material. As for CLVs, for any time $t$ , $\boldsymbol{v}^{\bot }(t)$ lives in the same function space as the primal solution $\boldsymbol{u}(t)$ , hence $\boldsymbol{v}^{\bot }$ also has $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D70C}\boldsymbol{U}$ and $\unicode[STIX]{x1D70C}E$ components. In figure 12, for each parameter, we plot the $|\unicode[STIX]{x1D70C}\boldsymbol{U}|$ field of $\boldsymbol{v}^{\bot }(t)$ . Unlike CLVs, both the direction and the magnitude of shadowing directions are meaningful, since shadowing directions reflect not only in what direction flow fields will change due to perturbations in system parameters $s$ , but also how large the change will be. Hence, we normalize shadowing directions by free-stream constants and unit parameters to preserve information about magnitudes.

Figure 12. Magnitude of the $\unicode[STIX]{x1D70C}\boldsymbol{U}$ component of the shadowing direction $\boldsymbol{v}^{\bot }=\text{d}\boldsymbol{u}^{\bot }/\text{d}s$ on the finer mesh, which is computed by FD-NILSS on a trajectory of time length $T=158t_{0}$ : (a)  $U$ as the parameter, normalized by $\unicode[STIX]{x1D70C}_{0}U_{0}/U_{0}=\unicode[STIX]{x1D70C}_{0}$ ; (b)  $\unicode[STIX]{x1D714}$ as the parameter, normalized by $\unicode[STIX]{x1D70C}_{0}U_{0}/\unicode[STIX]{x1D714}_{0}$ ; and (c) primal flow field, normalized by $\unicode[STIX]{x1D70C}_{0}U_{0}$ . Left: snapshots at $0.5T$ ; right: averaged from $66t_{0}$ to $93t_{0}$ . Notice that the left of (a) and (b) are coloured by log scale.

If we perturb system parameters by $\unicode[STIX]{x0394}U$ , the shadowing solution is a flow field with free-stream property $\unicode[STIX]{x1D70C}_{0}(U_{0}+\unicode[STIX]{x0394}U)$ . Since $\boldsymbol{v}^{\bot }$ points from the base solution to the shadowing solution, we expect the free-stream area of $\boldsymbol{v}^{\bot }$ to be $\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x0394}U/\unicode[STIX]{x0394}U=\unicode[STIX]{x1D70C}_{0}$ . The above intuition is verified on the left of figure 12(a): with $U$ as the system parameter, the free-stream area of the $\unicode[STIX]{x1D70C}\boldsymbol{U}$ component of the shadowing direction $\boldsymbol{v}^{\bot }$ has magnitude $1$ , since we are normalizing by $\unicode[STIX]{x1D70C}_{0}$ . On the other hand, perturbing rotation speed $\unicode[STIX]{x1D714}$ should have little effect on the free-stream area: this is also verified in figure 12(b), where $\boldsymbol{v}^{\bot }$ has small magnitude in the free-stream area.

We observe by comparing the left of figure 12(a) and figure 12(b) that, in the wake area, both $\boldsymbol{v}^{\bot }$ have similar magnitude. This fact is further confirmed by the fact that in figure 11 both $\Vert \boldsymbol{v}^{\bot }\Vert$ are similarly large after normalization. This fits our intuition that, if we perturb cylinder rotation by $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}=0.01\unicode[STIX]{x1D714}_{0}$ , then the circumferential speed is changed by $0.01\unicode[STIX]{x1D714}_{0}D/2=0.01U_{0}/2$ : this extra rotation should have similar magnitude of impact on the wake as changing the free-stream velocity by $\unicode[STIX]{x0394}U=0.01U_{0}$ .

We compare the shadowing directions with the primal flow field in figure 12. On one hand, shadowing directions are correlated to the primal flow field, since the peak values of shadowing directions are typically located close to the vortices, indicating that vortices structures are the most sensitive to parameter perturbations.

On the other hand, because shadowing directions are perturbations, they are very different from the primal flow. The primal flow has several structures that are insensitive to parameter perturbations, whereas parameter perturbations may bring new structures to the primal flow. For example, when $\unicode[STIX]{x1D714}$ is parameter, the stagnant area in front of the cylinder in primal flow changes to a smooth halo around the cylinder in the shadowing directions, due to the fact that perturbing $\unicode[STIX]{x1D714}$ leads to a uniform acceleration around the cylinder. Moreover, the shadowing direction has a larger active area than the primal flow. An extreme example is the left of figure 12(a), where the entire free stream of shadowing directions has value one, as we discussed. Indeed, the convolution formula of shadowing directions shows that shadowing directions are affected by both stable and unstable CLVs (Ni Reference Ni2018b ), the union of whose active areas cover the entire flow field.

4.5 Results of sensitivities

With shadowing directions, we can compute sensitivities of long-time-averaged objectives from (4.13) and (4.14) with little extra cost. In fact, as we shall see, using only a trajectory of time length 106 $t_{0}$ should be enough to compute accurate sensitivities. In this subsection, we investigate the effect of $U$ on the averaged drag force $\langle D_{r}\rangle$ . We will normalize $\langle D_{r}\rangle$ by $F_{0}=0.5\unicode[STIX]{x1D70C}U_{0}^{2}DZ=8.031\times 10^{-5}$ . The objective for system parameter $\unicode[STIX]{x1D714}$ is averaged lift $\langle L\rangle$ , which is also normalized by  $F_{0}$ .

We plot the history of the sensitivities computed via shadowing directions in figure 13. For each $T$ , we run NILSS from the same initial condition of primal and tangent solutions to $T$ , solve the shadowing direction, using which we compute the sensitivity corresponding to $T$ . The uncertainty of the final sensitivity is estimated by the smallest interval which bounds the history of the sensitivity and whose size shrinks as $T^{-0.5}$ . We can see that for both meshes the uncertainty of sensitivities is below 10 %.

Figure 13. History plots of sensitivities computed by FD-NILSS: (a) sensitivity of the drag with respect to $U$ , and (b) sensitivity of the lift with respect to $\unicode[STIX]{x1D714}$ .

Figure 14. Sensitivities for different choices of parameters and objectives: (a) drag versus $U$ , and (b) lift versus $\unicode[STIX]{x1D714}$ . Vertical bars indicate 95 % confidence intervals of the averaged objectives. The straight lines though the middle data points indicate the sensitivities computed by FD-NILSS.

We explain how to compute the uncertainty of long-time-averaged objectives, which is approximated by $\langle J\rangle _{T^{\prime }}$ , the averaged instantaneous objective over $T^{\prime }=7.68\times 10^{-3}=1014t_{0}$ . To get the uncertainty due to taking the average over finite time, we divide the history of $J(t)$ into five equally long parts. We denote the objectives averaged over each of the five parts by $J_{1},\ldots ,J_{5}$ , whose corrected sample standard deviation is denoted by $\unicode[STIX]{x1D70E}^{\prime }$ . We assume that the standard deviation of $\langle J\rangle _{T^{\prime }}$ is proportional to ${T^{\prime }}^{-0.5}$ , so we use $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}^{\prime }/\sqrt{5}$ as the standard deviation of $\langle J\rangle _{T^{\prime }}$ . We further assume $\pm 2\unicode[STIX]{x1D70E}$ yields the 95 % confidence interval, which is indicated by the vertical bar in figure 14.

The sensitivities computed via shadowing directions are shown in figure 14. For both meshes our sensitivities correctly reflect the trend between averaged objectives and system parameters. This confirms the physical meaning of the shadowing direction: it reveals the sensitivity of all physical phenomena in this flow field with respect to perturbations in parameters, in such a way that the average of the sensitivity is the sensitivity of the average.

The cost of FD-NILSS for computing a sensitivity is mainly in integrating the primal solution over $400\times 200\times 32=2.6\times 10^{6}$ time steps. Here 400 is the number of segments, 200 is the number of time steps in each segment, and 32 is the number of primal solutions computed. In FD-NILSS we need one inhomogeneous tangent and 30 homogeneous tangent solutions, and each tangent solution is approximated by the difference between a perturbed solution and the same base solution: that is 32 primal solutions in total. On the other hand, a sensitivity can be revealed by some regression among the five pairs of objectives and parameters in figure 14, the total cost of which is $3.8\times 10^{6}$ steps of primal simulation: this cost is similar to that of FD-NILSS.

NILSS can be further accelerated when we have a tangent solver, by taking advantage of the fact that all adjoint solutions use the same Jacobian $\unicode[STIX]{x2202}_{\boldsymbol{u}}\,\boldsymbol{f}$ . That is, we can integrate all tangent solutions simultaneously without repeatedly loading $\unicode[STIX]{x2202}_{\boldsymbol{u}}\,\boldsymbol{f}$ into the computer CPU, which is the most time-consuming procedure in the numerical integration. In other words, we can perform, at each time step, one matrix–matrix product, where the second matrix is composed of several tangent solutions, rather than performing several matrix–vector products, where we need to reload the matrix for each product. (This idea was brought up during private discussions with Pablo Fernandez, who co-authored FD-NILSS.)

Additionally, as discussed in Ni & Wang (Reference Ni and Wang2017) and Ni et al. (Reference Ni, Wang, Fernandez and Talnikar2018), FD-NILSS and NILSS have smaller marginal cost for new parameters; for cases with more parameters than objectives or cases when we have only adjoint solvers, the non-intrusive least-squares adjoint shadowing (NILSAS) method (Ni Reference Ni2018a ) can compute the gradient of one objective with respect to many parameters in one run.

From the above discussions, together with the implication of the conjecture in § 3.5 that there should be a small fraction of unstable CLVs for open flows, we believe that NILSS and NILSAS are competitive in terms of efficiency among all sensitivity algorithms.

As we can see, although our system is not uniform hyperbolic, the sensitivities computed by shadowing methods are still valid for both meshes. We thus add one more piece of evidence to the chaotic hypothesis: that is, although shadowing methods are logically consequential to uniform hyperbolicity, they are still valid for our fluid system, which only exhibits some hyperbolic phenomena (angles between far-apart CLVs are large), but is not uniform hyperbolic.