2.1 Lyapunov exponents

The thermoacoustic problem is governed by partial differential equations, i.e. the compressible Navier–Stokes equation with equations for the chemistry, and mass and energy conservation. After spatial discretisation, the thermoacoustic problem is formally an autonomous dynamical system

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\dot{\boldsymbol{q}}(t)=\boldsymbol{F}(\boldsymbol{q}(t)),\\ \boldsymbol{q}(0)=\boldsymbol{q}_{0},\end{array}\right\}\end{eqnarray}$$

where the overdot $\dot{(\,)}$ is Newton’s notation for time differentiation; $\boldsymbol{q}\in \mathbb{R}^{N}$ is the state vector (e.g. pressure and velocity at each discrete location), where the integer $N$ denotes the discrete degrees of freedom; the subscript $0$ denotes the initial condition; and $\boldsymbol{F}:\mathbb{R}^{N}\rightarrow \mathbb{R}^{N}$ is a nonlinear smooth function, which encapsulates the discretised boundary conditions. We are interested in the evolution of small perturbations, therefore we split the solution as

(2.2) $$\begin{eqnarray}\displaystyle \boldsymbol{q}(t)=\bar{\boldsymbol{q}}(t)+\boldsymbol{q}^{\prime }(t), & & \displaystyle\end{eqnarray}$$

where $\bar{\boldsymbol{q}}(t)$ is the unperturbed solution of (2.1) such that $\Vert \bar{\boldsymbol{q}}(t)\Vert \sim O(1)$ , and $\boldsymbol{q}^{\prime }(t)$ is the small perturbation such that $\Vert \bar{\boldsymbol{q}}^{\prime }(t)\Vert \sim O(\unicode[STIX]{x1D716})$ , where $\unicode[STIX]{x1D716}\rightarrow 0$ . The perturbation is governed by the tangent equation

(2.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\dot{\boldsymbol{q}}^{\prime }=\unicode[STIX]{x1D645}(t)\boldsymbol{q}^{\prime },\\ \boldsymbol{q}^{\prime }(0)=\boldsymbol{q}_{0}^{\prime },\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D645}(t)\equiv \text{d}\boldsymbol{F}/\text{d}\boldsymbol{q}|_{\bar{\boldsymbol{q}}(t)}$ is the Jacobian. To define the Lyapunov exponents, it is convenient to introduce the tangent propagator, which maps the perturbation, $\boldsymbol{q}^{\prime }$ , from time $t$ to time $\tilde{t}$ , as

(2.4) $$\begin{eqnarray}\displaystyle \boldsymbol{q}^{\prime }(t+\tilde{t})=\unicode[STIX]{x1D648}(t,\tilde{t})\boldsymbol{q}^{\prime }(t). & & \displaystyle\end{eqnarray}$$

The tangent propagator is governed by the matrix equation

(2.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\text{d}\unicode[STIX]{x1D648}}{\text{d}\tilde{t}}=\unicode[STIX]{x1D645}(\tilde{t})\unicode[STIX]{x1D648},\\ \unicode[STIX]{x1D648}(t,0)=\unicode[STIX]{x1D644},\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity matrix. Setting $t=0$ without loss of generality, the norm of an infinitesimal perturbation, $\boldsymbol{q}_{0}^{\prime }$ , to the initial condition, $\bar{\boldsymbol{q}}_{0}$ , asymptotically grows (or decays) as

(2.6) $$\begin{eqnarray}\Vert \boldsymbol{q}^{\prime }(\tilde{t})\Vert \cong \Vert \boldsymbol{q}_{0}^{\prime }\Vert \text{e}^{\unicode[STIX]{x1D706}(\boldsymbol{q}_{0}^{\prime },\bar{\boldsymbol{q}}_{0})\tilde{t}},\end{eqnarray}$$

where $\cong$ means ‘asymptotically equal to’, and

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D706}(\boldsymbol{q}_{0}^{\prime },\bar{\boldsymbol{q}})=\lim _{\tilde{t}\rightarrow \infty }\frac{1}{\tilde{t}}\log \frac{\Vert \unicode[STIX]{x1D648}(0,\tilde{t})\boldsymbol{q}_{0}^{\prime }\Vert }{\Vert \boldsymbol{q}_{0}^{\prime }\Vert }\end{eqnarray}$$

is the characteristic Lyapunov exponent. Oseledets’ theorem (Oseledets Reference Oseledets1968) shows that there exist $m\leqslant N$ distinct Lyapunov exponents $\unicode[STIX]{x1D706}_{1}(\bar{\boldsymbol{q}})>\unicode[STIX]{x1D706}_{2}(\bar{\boldsymbol{q}})>\cdots >\unicode[STIX]{x1D706}_{m}(\bar{\boldsymbol{q}})$ , which provide a filtration of the tangent space ${\mathcal{T}}_{\bar{\boldsymbol{q}}}$ , into subspaces $S_{i}$ , i.e. ${\mathcal{T}}_{\bar{\boldsymbol{q}}}\equiv S_{1}\supset S_{2}\supset \cdots \supset S_{m}$ , such that $\boldsymbol{q}_{0}^{\prime }\in S_{j}\backslash S_{j+1}\Leftrightarrow \unicode[STIX]{x1D706}(\boldsymbol{q}_{0}^{\prime },\bar{\boldsymbol{q}})=\unicode[STIX]{x1D706}_{j}(\bar{\boldsymbol{q}})$ . Furthermore, Oseledets’ theorem shows that $\unicode[STIX]{x1D706}_{j}(\bar{\boldsymbol{q}})$ are constants of the attractor $\bar{\boldsymbol{q}}$ , and, in ergodic systems, $\unicode[STIX]{x1D706}_{j}$ do not depend on the initial condition, $\bar{\boldsymbol{q}}_{0}$ . Physically, the Lyapunov exponents are the average exponential contraction/expansion rates of an infinitesimal volume of the phase space moving along the attractor. For example, as shown in § 3.1, if the attractor is a fixed point, the Lyapunov exponents are equal to the real part of the eigenvalues of the Jacobian at the fixed point. Similarly, as shown in § 3.2, if the attractor is a limit cycle, the Lyapunov exponents are equal to the real part of the Floquet exponents.

2.2 Oseledets’ splitting and covariant Lyapunov vectors

The Lyapunov exponents are invariant measures of the attractor, however, they do not inform on the directions along which the infinitesimal volume of the phase space contracts/expands. Such directions are provided by the Oseledets splitting, which is composed by the Lyapunov subspaces, which are, in turn, spanned by the covariant Lyapunov vectors. First, the Oseledets matrix (Oseledets Reference Oseledets1968) is defined as

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D729}^{\pm }(t)=\lim _{\tilde{t}\rightarrow \pm \infty }\frac{1}{2\tilde{t}}\log [\unicode[STIX]{x1D648}(t,\tilde{t})^{\text{T}}\unicode[STIX]{x1D648}(t,\tilde{t})].\end{eqnarray}$$

This matrix is called ‘forward’ if $\tilde{t}\rightarrow +\infty$ or ‘backward’ if $\tilde{t}\rightarrow -\infty$ . The spectrum of this matrix contains $m\leqslant N$ distinct eigenvalues, which are the Lyapunov exponents of the system $\unicode[STIX]{x1D706}_{1}>\unicode[STIX]{x1D706}_{2}>\cdots >\unicode[STIX]{x1D706}_{m}$ . However, the eigenvectors of the forward and backward matrices differ from each other and are not invariant under time reversal. To gain more insight in the Oseledets matrix, consider the singular value decomposition $\unicode[STIX]{x1D648}(t,\tilde{t})=\unicode[STIX]{x1D650}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D651}^{\text{T}}$ , where $\unicode[STIX]{x1D650}$ and $\unicode[STIX]{x1D651}$ are orthogonal matrices and $\unicode[STIX]{x1D72E}$ is a diagonal matrix with non-negative real entries (the singular values). We can obtain an eigenvalue decomposition of the argument of the logarithm of (2.8), $\unicode[STIX]{x1D648}^{\text{T}}\unicode[STIX]{x1D648}=\unicode[STIX]{x1D651}(\unicode[STIX]{x1D72E}^{\text{T}}\unicode[STIX]{x1D72E})\unicode[STIX]{x1D651}^{\text{T}}=\unicode[STIX]{x1D651}\unicode[STIX]{x1D72E}^{2}\unicode[STIX]{x1D651}^{\text{T}}$ , which, after applying the logarithm, becomes $\unicode[STIX]{x1D651}\log (\unicode[STIX]{x1D72E}^{2})\unicode[STIX]{x1D651}^{\text{T}}=2\unicode[STIX]{x1D651}\log (\unicode[STIX]{x1D72E})\unicode[STIX]{x1D651}^{\text{T}}$ . Thus, equation (2.8) can be rewritten as $\unicode[STIX]{x1D729}^{\pm }(t)=\lim _{\tilde{t}\rightarrow \pm \infty }\unicode[STIX]{x1D651}(\log (\unicode[STIX]{x1D72E}(t,\tilde{t}))/\tilde{t})\unicode[STIX]{x1D651}^{\text{T}}$ , which shows that the eigenvalues of $\unicode[STIX]{x1D729}^{\pm }$ are the Lyapunov exponents, which are equal to the exponential average of the singular values of $\unicode[STIX]{x1D648}(t,\tilde{t})$ . Let $V_{j}^{\pm }(t)$ be the $j$ th eigenspace of the forward (backward) Oseledets matrix, then the Oseledets subspaces are defined as

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}_{j}^{+}(t)=V_{j}^{+}(t)\oplus \cdots \oplus V_{m}^{+}(t), & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}_{j}^{-}(t)=V_{1}^{-}(t)\oplus \cdots \oplus V_{j}^{-}(t), & \displaystyle\end{eqnarray}$$

where $\oplus$ is the direct sum. The Oseledets subspaces have the property

(2.11) $$\begin{eqnarray}\lim _{\tilde{t}\rightarrow \infty }\frac{1}{\tilde{t}}\log \frac{\Vert \unicode[STIX]{x1D648}(t,\tilde{t})\boldsymbol{q}^{\prime }(t)\Vert }{\Vert \boldsymbol{q}^{\prime }(t)\Vert }=\unicode[STIX]{x1D706}_{j},\quad \text{for }\boldsymbol{q}^{\prime }(t)\in \unicode[STIX]{x1D6E4}_{j}^{\pm }(t)\backslash \unicode[STIX]{x1D6E4}_{j+1}^{\pm }(t)\end{eqnarray}$$

and a nested structure $\mathbb{R}^{N}=\unicode[STIX]{x1D6E4}_{1}^{+}(t)\supset \unicode[STIX]{x1D6E4}_{2}^{+}(t)\supset \cdots \supset \unicode[STIX]{x1D6E4}_{m}^{+}(t)\supset \unicode[STIX]{x1D6E4}_{m+1}^{+}(t)\equiv \emptyset$ and $\mathbb{R}^{N}=\unicode[STIX]{x1D6E4}_{m}^{-}(t)\supset \unicode[STIX]{x1D6E4}_{m-1}^{-}(t)\supset \cdots \supset \unicode[STIX]{x1D6E4}_{1}^{-}(t)\supset \unicode[STIX]{x1D6E4}_{0}^{-}(t)\equiv \emptyset$ . By intersecting the Oseledets subspaces, we obtain the Lyapunov subspaces

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{j}(t)=\unicode[STIX]{x1D6E4}_{j}^{+}(t)\cap \unicode[STIX]{x1D6E4}_{j}^{-}(t),\quad j\in \{1,\ldots ,m\},\end{eqnarray}$$

which compose the Oseledets splitting. The Lyapunov subspaces are (i) generally non-orthogonal to each other, (ii) covariant with the dynamics, i.e. $\unicode[STIX]{x1D648}(t,\tilde{t})\unicode[STIX]{x1D6FA}_{j}(t)=\unicode[STIX]{x1D6FA}_{j}(t+\tilde{t})$ , and (iii) invariant under time reversal. Each vector $\unicode[STIX]{x1D753}_{j}(t)$ of a set that spans one of the Lyapunov subspaces is a covariant Lyapunov vector associated with the Lyapunov exponent $\unicode[STIX]{x1D706}_{j}$ . If a trajectory is infinitesimally perturbed at some time $t_{1}$ along a covariant Lyapunov vector, the perturbation will grow at an exponential rate dictated by the associated Lyapunov exponent and will stay aligned with that same covariant Lyapunov vector (figure 1).

Figure 1. Schematic diagram of covariant Lyapunov vectors and perturbations on an unperturbed trajectory (solid grey line). Three covariant Lyapunov vectors are shown at two different instants, each associated with a different Lyapunov exponent, which can be positive, zero or negative. The decay/growth of three perturbations, along the stable (green), neutral (orange), unstable (red) covariant Lyapunov vector, respectively. The resulting perturbed trajectories (dashed lines), converge, remain at a constant distance, or diverge, respectively, to/from the unperturbed trajectory, depending on the sign of the Lyapunov exponent. This explains why trajectories emanating from two very close initial conditions will almost surely diverge in chaotic systems – it is highly unlikely for the vector connecting the two initial conditions not to have a component in the direction of the unstable covariant Lyapunov vector.

We derive the equation that governs the covariant Lyapunov vectors. First, because the Lyapunov subspaces are covariant with the dynamics, the following definition holds

(2.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D648}(t,\tilde{t})\unicode[STIX]{x1D753}(t)=\unicode[STIX]{x1D702}(t,\tilde{t})\unicode[STIX]{x1D753}(t+\tilde{t}), & & \displaystyle\end{eqnarray}$$

where

(2.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D648}(t,\tilde{t})=\text{e}^{\int _{t}^{t+\tilde{t}}\unicode[STIX]{x1D645}(\unicode[STIX]{x1D712})\,\text{d}\unicode[STIX]{x1D712}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D702}(t,\tilde{t})$ is a scalar that measures the asymptotic growth of the norm and allows $\unicode[STIX]{x1D753}(t+\tilde{t})$ to have any desired bounded norm. Substituting (2.14) in (2.13) and differentiating with respect to $\tilde{t}$ results in

(2.15) $$\begin{eqnarray}\unicode[STIX]{x1D645}(t+\tilde{t})\text{e}^{\int _{t}^{t+\tilde{t}}\unicode[STIX]{x1D645}(\unicode[STIX]{x1D712})\,\text{d}\unicode[STIX]{x1D712}}\unicode[STIX]{x1D753}(t)=\frac{\text{d}\unicode[STIX]{x1D702}(t,\tilde{t})}{\text{d}\tilde{t}}\unicode[STIX]{x1D753}(t+\tilde{t})+\unicode[STIX]{x1D702}(t,\tilde{t})\frac{\text{d}\unicode[STIX]{x1D753}(t+\tilde{t})}{\text{d}\tilde{t}}.\end{eqnarray}$$

By setting $t=0$ and omitting the explicit dependence on $t=0$ , we obtain

(2.16) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D753}}{\text{d}\tilde{t}}=\unicode[STIX]{x1D645}(\tilde{t})\unicode[STIX]{x1D753}(\tilde{t})-\frac{1}{\unicode[STIX]{x1D702}(\tilde{t})}\frac{\text{d}\unicode[STIX]{x1D702}(\tilde{t})}{\text{d}\tilde{t}}\unicode[STIX]{x1D753}(\tilde{t}),\end{eqnarray}$$

for any $\unicode[STIX]{x1D702}(\tilde{t})\not =0$ . Moreover, we know from Oseledets’ theorem that

(2.17) $$\begin{eqnarray}\Vert \unicode[STIX]{x1D648}(\tilde{t})\unicode[STIX]{x1D753}(0)\Vert \cong \text{e}^{\unicode[STIX]{x1D706}\tilde{t}}\Vert \unicode[STIX]{x1D753}(0)\Vert ,\end{eqnarray}$$

which shows that $\unicode[STIX]{x1D702}(\tilde{t})\cong \text{e}^{\unicode[STIX]{x1D706}\tilde{t}}$ . If we choose to have a bounded non-zero covariant Lyapunov vector, i.e. $0<\Vert \unicode[STIX]{x1D753}\Vert <\infty$ , equation (2.16) becomes

(2.18) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D753}}{\text{d}\tilde{t}}=\unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D753}.\end{eqnarray}$$

It is easier to mathematically manipulate and numerically solve (2.18) than (2.12). Moreover, equation (2.18) provides a clear picture of the evolution of a covariant Lyapunov vector: the vector is evolved by the tangent dynamics $\unicode[STIX]{x1D645}\unicode[STIX]{x1D753}$ , while the extra term $-\unicode[STIX]{x1D706}\unicode[STIX]{x1D753}$ guarantees that its norm is bounded. It can be shown that if the attractor is periodic or chaotic, there is a neutral mode ( $\unicode[STIX]{x1D706}=0$ ), where $\unicode[STIX]{x1D753}$ is collinear to $\dot{\bar{\boldsymbol{q}}}$ (Katok & Hasselblatt Reference Katok and Hasselblatt1995).

In the remainder of this paper, $t=0$ without loss of generality and the tilde, $\tilde{()}$ , is dropped for brevity.

2.3 Numerical computation of Lyapunov exponents and covariant Lyapunov vectors

We use a robust algorithm (Ginelli et al. Reference Ginelli, Poggi, Turchi, Chaté, Livi and Politi2007, Reference Ginelli, Chaté, Livi and Politi2013), called the QR algorithm for brevity, to calculate the Lyapunov spectrum and covariant vectors. The algorithm evolves a set of $m$ column vectors $\boldsymbol{g}_{j},\,j=\{1,\ldots ,m\}$ of an $N\times m$ matrix $\unicode[STIX]{x1D648}$ , via the tangent equation (2.5). Because the $\boldsymbol{g}_{j}$ will likely have a component in the most unstable (least stable) space $\unicode[STIX]{x1D6FA}_{1}$ , their norm will exponentially grow (decay) at rate $\unicode[STIX]{x1D706}_{1}$ , which is bound to numerically overflow (underflow). To overcome this numerical instability, the QR algorithm executes periodic orthonormalisations of $\unicode[STIX]{x1D648}$ . By denoting the time step with a superscript, the calculation of the Lyapunov exponents and covariant Lyapunov vectors is enabled by the following algorithm.

1. (i) Set the initial condition $\bar{\boldsymbol{q}}^{0}$ and initialise $\unicode[STIX]{x1D648}^{0}$ to a random orthonormal set of vectors $[\boldsymbol{g}_{1}\ldots \boldsymbol{g}_{m}]$ .

2. (ii) Evolve $\bar{\boldsymbol{q}}^{0},\unicode[STIX]{x1D648}^{0}$ using (2.1), (2.5) for $n_{su}$ iterations, where $n_{su}$ is called the spinup time, which must be sufficiently large such that $\bar{\boldsymbol{q}}^{n_{su}}$ is in the attractor (to some numerical tolerance).

3. (iii) Evolve $\bar{\boldsymbol{q}}^{j},\unicode[STIX]{x1D648}^{j}$ for $n_{QR}$ iterations.

4. (iv) Perform QR decomposition on $\unicode[STIX]{x1D648}^{j}$ , obtaining $\unicode[STIX]{x1D64C}^{j},\boldsymbol{R}^{j}$ . Store $\unicode[STIX]{x1D64C}^{j},\boldsymbol{R}^{j}$ , and set $\unicode[STIX]{x1D648}^{j}:=\unicode[STIX]{x1D64C}^{j}$ . If $j<n_{T}$ , where $n_{T}$ is the total number of iterations corresponding to the total simulation time, go back to Item (iii).

5. (v) Randomly initialise an upper triangular matrix $\unicode[STIX]{x1D63E}^{n_{T}}$ of the same dimension as all $\boldsymbol{R}^{j}$ .

6. (vi) Evolve $\unicode[STIX]{x1D63E}^{j}$ backward by solving $\boldsymbol{R}^{j}\unicode[STIX]{x1D63E}^{j}=\unicode[STIX]{x1D63E}^{j+1}$ for $\unicode[STIX]{x1D63E}^{j}$ and subsequently normalise its columns, i.e. ensuring $\sum _{k}(\unicode[STIX]{x1D60A}_{lk}^{j})^{2}=1$ .

7. (vii) Compute Lyapunov exponents: $[\unicode[STIX]{x1D706}_{1}\ldots \unicode[STIX]{x1D706}_{m}]=((n_{T}-n_{su})\unicode[STIX]{x0394}t)^{-1}\sum _{j=n_{su}}^{n_{T}}\log (|\text{diag}(\boldsymbol{R}^{j})|)$ , where $\unicode[STIX]{x0394}t$ is the time step.

8. (viii) Compute covariant Lyapunov vectors: $[\unicode[STIX]{x1D753}_{1}|\ldots |\unicode[STIX]{x1D753}_{m}]^{j}=\unicode[STIX]{x1D64C}^{j}\unicode[STIX]{x1D63E}^{j}$ , valid only for $j\in [n_{su},n_{T}-n_{sd}]$ , where $n_{sd}$ is the spindown time, which must be sufficiently large for $\unicode[STIX]{x1D63E}^{j}$ to converge to the covariant Lyapunov vector expansion coefficients.

2.4 Hyperbolicity

A strange attractor is hyperbolic if there is a splitting of the tangent space into stable, neutral and unstable subspaces at every point of the trajectory, $\bar{\boldsymbol{q}}(t)$ . Formally, ${\mathcal{T}}_{\bar{\boldsymbol{q}}}=E_{\bar{\boldsymbol{q}}}^{s}\oplus E_{\bar{\boldsymbol{q}}}^{n}\oplus E_{\bar{\boldsymbol{q}}}^{u}$ , where $E_{\bar{\boldsymbol{q}}}^{s}$ and $E_{\bar{\boldsymbol{q}}}^{u}$ are the stable and unstable subspaces of dimension $N^{s}$ and $N^{u}$ , defined by the directions along which the derivative contracts and expands, respectively, and $E_{\bar{\boldsymbol{q}}}^{n}$ is the one-dimensional neutral subspace. (Consequently, a quasi-periodic solution is not hyperbolic because it has at least two zero Lyapunov exponents, i.e. $E_{\bar{\boldsymbol{q}}}^{n}$ is at least two-dimensional.) Hyperbolicity has profound implications on the behaviour of a dynamical system. The existence of unstable subspaces gives rise to exponentially diverging trajectories, which in turn gives rise to unpredictable dynamics in the long term. Furthermore, hyperbolicity typically implies structural stability of the attractor, i.e. the qualitative behaviour of the attractor does not change if the system is slightly perturbed. In the problem of computing sensitivities, hyperbolicity is crucial because it determines whether the time-averaged cost functional responds smoothly to perturbations to the parameters (Ruelle Reference Ruelle1980). Indeed, the most robust sensitivity algorithms (e.g. Wang Reference Wang2013; Blonigan & Wang Reference Blonigan and Wang2014; Wang et al. Reference Wang, Hu and Blonigan2014; Blonigan Reference Blonigan2017; Ni & Wang Reference Ni and Wang2017) rely on the shadowing lemma (Bowen & Ruelle Reference Bowen and Ruelle1975; Pilyugin Reference Pilyugin2006), which is valid only in hyperbolic systems. Importantly, it has been hypothesised by Gallavotti & Cohen (Reference Gallavotti and Cohen1995), Gallavotti (Reference Gallavotti2006) that most physical dynamical systems develop asymptotically on an attracting set, the dynamics of which can be regarded as hyperbolic. This is called the chaotic hypothesis, which stems from the measure theory of turbulence of Ruelle (Reference Ruelle1980). In order to verify hyperbolicity in a numerical simulation, the method described in Takeuchi et al. (Reference Takeuchi, Yang, Ginelli, Radons and Chaté2011) is used here. The angles between the three pairs of subspaces, $\unicode[STIX]{x1D703}_{u,n}=\angle (E_{\bar{\boldsymbol{q}}}^{u},E_{\bar{\boldsymbol{q}}}^{n})$ , $\unicode[STIX]{x1D703}_{u,s}=\angle (E_{\bar{\boldsymbol{q}}}^{u},E_{\bar{\boldsymbol{q}}}^{s})$ , $\unicode[STIX]{x1D703}_{n,s}=\angle (E_{\bar{\boldsymbol{q}}}^{n},E_{\bar{\boldsymbol{q}}}^{s})$ , are computed. These angles are computed by using the principal angles, $\cos (\unicode[STIX]{x1D703}_{A,B})=\bar{\unicode[STIX]{x1D70E}}(\unicode[STIX]{x1D64C}_{A}\unicode[STIX]{x1D64C}_{B})$ , where matrices $\unicode[STIX]{x1D64C}_{A}$ and $\unicode[STIX]{x1D64C}_{B}$ define the orthonormal bases of any subspaces $A$ and $B$ (not only $E_{\bar{\boldsymbol{q}}}^{u}$ , $E_{\bar{\boldsymbol{q}}}^{n}$ , $E_{\bar{\boldsymbol{q}}}^{s}$ ), respectively, and $\bar{\unicode[STIX]{x1D70E}}$ is the largest singular value. Then, a probability density function of each of these angles is extracted via a histogram of the time series (e.g. figure 15 b). The system behaves hyperbolically if there are no tangencies between the subspaces, i.e. the value of the probability density functions at $\unicode[STIX]{x1D703}=0$ is 0.

2.5 Shadowing lemma

Shadowing-based sensitivity methods are centred around the shadowing lemma. The shadowing lemma exists both for discrete or continuous dynamical systems, but we will only present its discrete version because most engineering problems are numerically discretised to be solved.

Definition 1 ( $\unicode[STIX]{x1D716}$ -pseudo-orbit).

An $\unicode[STIX]{x1D716}$ -pseudo-orbit for the map $\boldsymbol{f}$ is a sequence of points $\{\boldsymbol{y}_{n}\}$ , such that

(2.19) $$\begin{eqnarray}\Vert \boldsymbol{y}_{n+1}-\boldsymbol{f}(\,\boldsymbol{y}_{n})\Vert <\unicode[STIX]{x1D716}.\end{eqnarray}$$

An $\unicode[STIX]{x1D716}$ -pseudo-orbit is thus a series $\{\boldsymbol{y}_{n}\}$ where each point $\boldsymbol{y}_{n+1}$ is at most $\unicode[STIX]{x1D716}$ away from the true iterate $\boldsymbol{f}(\,\boldsymbol{y}_{n})$ of the previous point $\boldsymbol{y}_{n}$ .

Definition 2 ( $\unicode[STIX]{x1D6FF}$ -shadow-orbit).

An actual orbit $\{\boldsymbol{x}_{n}\}$ , where $\boldsymbol{x}_{n}=\boldsymbol{f}^{n}(\boldsymbol{x}_{0})$ , is said to be a $\unicode[STIX]{x1D6FF}$ -shadow-orbit of the $\unicode[STIX]{x1D716}$ -pseudo-orbit $\{{\boldsymbol{y}_{n}\}}_{a<n<b}$ if

(2.20) $$\begin{eqnarray}\Vert \boldsymbol{x}_{n}-\boldsymbol{y}_{n}\Vert <\unicode[STIX]{x1D6FF}.\end{eqnarray}$$

In other words, an $\unicode[STIX]{x1D716}$ -pseudo-orbit is $\unicode[STIX]{x1D6FF}$ -shadowed by a true orbit if the true orbit is closer than $\unicode[STIX]{x1D6FF}$ from it at each point.

Intuitively, the shadowing lemma guarantees that, even in chaotic systems where infinitely close trajectories diverge, there is a trajectory of a slightly perturbed system that does not diverge from the unperturbed system.

2.6 Shadowing methods for sensitivity

The gradient of the time-averaged cost functional (1.1) explicitly reads

(2.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{\boldsymbol{s}}\langle {\mathcal{J}}\rangle \triangleq \frac{\text{d}}{\text{d}\boldsymbol{ s}}\left(\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{\text{T}}{\mathcal{J}}(\boldsymbol{s};\boldsymbol{q}(t))\,\text{d}t\right). & & \displaystyle\end{eqnarray}$$

In a chaotic attractor, the operations of differentiation and time average do not commute, i.e. $\unicode[STIX]{x1D735}_{\boldsymbol{s}}\langle {\mathcal{J}}\rangle \not =\langle \unicode[STIX]{x1D735}_{\boldsymbol{s}}{\mathcal{J}}\rangle$ , where

(2.22) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D735}_{\boldsymbol{s}}{\mathcal{J}}\rangle =\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{\text{T}}\left(\frac{\unicode[STIX]{x2202}{\mathcal{J}}}{\unicode[STIX]{x2202}\boldsymbol{s}}+\frac{\unicode[STIX]{x2202}{\mathcal{J}}}{\unicode[STIX]{x2202}\boldsymbol{q}}\frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}\boldsymbol{s}}\right)\,\text{d}t. & & \displaystyle\end{eqnarray}$$

Equation (2.22) is an unbounded quantity because $\unicode[STIX]{x2202}\boldsymbol{q}/\unicode[STIX]{x2202}\boldsymbol{s}$ lives in the tangent space, a subspace that is exponentially unstable in chaotic attractors. Shadowing-based sensitivity methods integrate the sensitivity of the cost functional along the shadow trajectory, which does not diverge from the attractor. This way, the quantity (2.22) is bounded and equal to $\unicode[STIX]{x1D735}_{\boldsymbol{s}}\langle {\mathcal{J}}\rangle$ . The original shadowing method (Wang Reference Wang2013) achieves this by calculating the perturbation to $\boldsymbol{F}$ in (2.1) due to a perturbation in a parameter $\boldsymbol{s}\rightarrow \boldsymbol{s}+\unicode[STIX]{x1D6FF}\boldsymbol{s}$ , that is, $\unicode[STIX]{x2202}\boldsymbol{F}/\unicode[STIX]{x2202}\boldsymbol{s}\cdot \unicode[STIX]{x1D6FF}\boldsymbol{s}$ . The perturbation is decomposed in the covariant Lyapunov vector basis of the baseline trajectory to obtain a set of independent ordinary differential equations, one for each mode. The solutions of these equations are the components of the shadow trajectory in the covariant Lyapunov vector basis. After obtaining the perturbed trajectory in the phase space, the sensitivities can be readily computed. A major drawback of the original shadowing method is the need to compute all the covariant Lyapunov vectors for all time steps, which is computationally expensive. The least-squares shadowing method (Wang et al. Reference Wang, Hu and Blonigan2014) overcomes this by finding a trajectory of the system at parameter value $\boldsymbol{s}+\unicode[STIX]{x1D6FF}\boldsymbol{s}$ that is close to a trajectory of the system at parameter $\boldsymbol{s}$ via solving a least squares minimisation problem, which minimises the distance between the two trajectories at regular checkpoints. While the least-squares shadowing method is faster than the original shadowing method, it still carries high computational cost, as it requires solving a linear system of dimension equal to the dimension of the phase space times the number of checkpoints. Ni & Wang (Reference Ni and Wang2017) developed the non-intrusive least-squares shadowing method, the computational cost of which scales only with the number of unstable covariant Lyapunov vectors. In this paper, we will apply the non-intrusive least-squares shadowing method to a chaotic thermoacoustic system.

3.1 Eigenvalue analysis of fixed points: connection with covariant Lyapunov vectors

Eigenvalue analysis determines the linear stability of a fixed point of $\boldsymbol{F}$ . Mathematically, in decomposition (2.2), $\bar{\boldsymbol{q}}$ does not depend on time. The linearised dynamics around the fixed point $\bar{\boldsymbol{q}}$ is governed by (2.3) where the Jacobian $\unicode[STIX]{x1D645}=\text{d}\boldsymbol{F}/\text{d}\boldsymbol{q}|_{\bar{\boldsymbol{q}}}$ is constant. The formal solution for an initial condition reads

(3.1) $$\begin{eqnarray}\boldsymbol{q}^{\prime }(t)=\text{e}^{\unicode[STIX]{x1D645}t}\boldsymbol{q}_{0}^{\prime }.\end{eqnarray}$$

By assuming that $\unicode[STIX]{x1D645}$ has a complete eigenbasis, i.e. it is not defective, $\boldsymbol{q}_{0}^{\prime }$ can be decomposed in the eigenbasis $\{\hat{\boldsymbol{q}}_{1},\ldots ,\hat{\boldsymbol{q}}_{N}\}$ , where $\hat{()}_{j}$ is an eigenvector of $\unicode[STIX]{x1D645}$ , as

(3.2) $$\begin{eqnarray}\boldsymbol{q}_{0}^{\prime }=\mathop{\sum }_{j=1}^{N}\underbrace{(\boldsymbol{q}_{0}^{\prime }\boldsymbol{\cdot }\hat{\boldsymbol{q}_{j}})}_{\triangleq \unicode[STIX]{x1D6FC}_{j}}\hat{\boldsymbol{q}}_{j},\end{eqnarray}$$

where, to keep a similar notation to covariant Lyapunov vector analysis, the eigenpairs are sorted in descending order according to the real part of the corresponding eigenvalue $\unicode[STIX]{x1D70E}_{j}$ , i.e. $j=1$ denotes the eigenpair with largest growth rate. Substituting (3.2) in (3.1) yields

(3.3) $$\begin{eqnarray}\displaystyle \boldsymbol{q}^{\prime }(t) & = & \displaystyle \text{e}^{\unicode[STIX]{x1D645}t}\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\hat{\boldsymbol{q}}_{j}\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\text{e}^{\unicode[STIX]{x1D70E}_{j}t}\hat{\boldsymbol{q}}_{j}.\end{eqnarray}$$

Substituting the perturbation (3.3) into the definition of Lyapunov exponent, equation (2.7), yields

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706} & = & \displaystyle \lim _{t\rightarrow \infty }\frac{1}{t}\log \frac{\displaystyle \left\Vert \mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\text{e}^{\unicode[STIX]{x1D70E}_{j}t}\hat{\boldsymbol{q}}_{j}\right\Vert }{\Vert \boldsymbol{q}_{0}^{\prime }\Vert }\nonumber\\ \displaystyle & = & \displaystyle \lim _{t\rightarrow \infty }\frac{1}{t}\log \left(\frac{|\text{e}^{\unicode[STIX]{x1D70E}_{k}t}|\,\Vert \unicode[STIX]{x1D6FC}_{k}\hat{\boldsymbol{q}}_{k}\Vert }{\Vert \boldsymbol{q}_{0}^{\prime }\Vert }\right)\nonumber\\ \displaystyle & = & \displaystyle \lim _{t\rightarrow \infty }\frac{1}{t}\left[\log (\text{e}^{\text{Re}(\unicode[STIX]{x1D70E}_{k})t})+\log \left(\frac{\Vert \unicode[STIX]{x1D6FC}_{k}\hat{\boldsymbol{q}}_{k}\Vert }{\Vert \boldsymbol{q}_{0}^{\prime }\Vert }\right)\right]\nonumber\\ \displaystyle & = & \displaystyle \text{Re}(\unicode[STIX]{x1D70E}_{k}),\end{eqnarray}$$

where $k$ is the first index such that $\unicode[STIX]{x1D6FC}_{k}\neq 0$ , and $\text{Re}$ denotes the real part. Equation (3.4) shows that the $k$ th Lyapunov exponent, $\unicode[STIX]{x1D706}_{k}$ , of a linear system on a fixed point is the real part of the eigenvalue of the Jacobian, $\text{Re}(\unicode[STIX]{x1D70E}_{k})$ . Physically, the Lyapunov exponent is the growth (or decay) rate of small perturbations on top of the steady solution. This means that the Lyapunov exponent associated with the perturbation $\boldsymbol{q}_{0}^{\prime }$ is $\text{Re}(\unicode[STIX]{x1D70E}_{k})$ if $\boldsymbol{q}_{0}^{\prime }$ does not belong to $\text{Span}(\hat{\boldsymbol{q}}_{1},\ldots ,\hat{\boldsymbol{q}}_{k-1})$ , which is the subspace spanned by the first $k-1$ eigenvectors. Because the covariant Lyapunov vector is an unsteady vector, we have to decompose it in a time-varying basis of the propagator $\text{e}^{\unicode[STIX]{x1D645}t}$ , which consists of the time-varying vectors $\hat{\boldsymbol{q}}_{j}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}$ (Curtain & Zwart Reference Curtain and Zwart1995), where $\text{i}$ is the imaginary unit. The angular frequency of the linear oscillation is denoted $\unicode[STIX]{x1D714}_{j}=\text{Im}(\unicode[STIX]{x1D70E}_{j})$ , where $\text{Im}$ is the imaginary part. Abusing notation by re-using $\unicode[STIX]{x1D6FC}_{j}$ as the coordinates in the new basis, we can decompose $\unicode[STIX]{x1D753}$ as

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D753}(t)=\mathop{\sum }_{j=1}^{N}\underbrace{(\unicode[STIX]{x1D753}_{0}\boldsymbol{\cdot }\hat{\boldsymbol{q}}_{j})}_{\triangleq \unicode[STIX]{x1D6FC}_{j}}\hat{\boldsymbol{q}}_{j}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t},\end{eqnarray}$$

such that its time derivative reads

(3.6) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D753}} & = & \displaystyle \mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}(\unicode[STIX]{x1D70E}_{j}-\text{Re}(\unicode[STIX]{x1D70E}_{j}))\hat{\boldsymbol{q}}_{j}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\unicode[STIX]{x1D645}\hat{\boldsymbol{q}}_{j}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}-\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D6FC}_{j}\hat{\boldsymbol{q}}_{j}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D6FC}_{j}\hat{\boldsymbol{q}}_{j}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}.\end{eqnarray}$$

In order to factor out $\unicode[STIX]{x1D706}_{j}$ , we consider the set of eigenvectors that share the same growth rate, $\text{Re}(\unicode[STIX]{x1D70E}_{j})$ , although they may have different angular frequencies, $\unicode[STIX]{x1D714}_{j}$ . Mathematically, by considering $\unicode[STIX]{x1D6FC}_{k}=0:k\notin \{\,j:\unicode[STIX]{x1D706}_{j}=\unicode[STIX]{x1D706}\}$ , the covariant Lyapunov vector equation (2.18) is recovered

(3.7) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D753}} & = & \displaystyle \unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D6FC}_{j}\hat{\boldsymbol{q}_{j}}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\unicode[STIX]{x1D706}\mathop{\sum }_{j,\unicode[STIX]{x1D6FC}_{j}\neq 0}\unicode[STIX]{x1D6FC}_{j}\hat{\boldsymbol{q}}_{j}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D753},\end{eqnarray}$$

which shows that covariant Lyapunov vectors and eigenvectors, grouped by growth rates, span the same subspaces. Furthermore, any real linear combination of $\hat{\boldsymbol{q}}_{j}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}$ is a covariant Lyapunov vector with the different $\hat{\boldsymbol{q}}_{j}$ corresponding to eigenvalues that have the same growth rate (i.e. same Lyapunov exponents). On the one hand, if $\unicode[STIX]{x1D70E}_{j}\in \mathbb{R}$ , there is only one such $j$ and thus $\unicode[STIX]{x1D753}=\hat{\boldsymbol{q}_{j}}$ is a covariant Lyapunov vector. (In principle, there may be cases where the spectrum contains one real eigenvalue and two complex conjugates with the same real part as the real eigenvalue. Although we do not consider this special case, the conclusions we draw still hold.) On the other hand, if $\unicode[STIX]{x1D70E}_{j}\in \mathbb{C}$ , there is a pair $\{\unicode[STIX]{x1D6FC}_{j}\hat{\boldsymbol{q}_{j}}\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t},\unicode[STIX]{x1D6FC}_{j}^{\ast }\hat{\boldsymbol{q}}_{j}^{\ast }\text{e}^{\text{i}\unicode[STIX]{x1D714}_{j}t}\}$ that defines a two-dimensional subspace in which all vectors are covariant Lyapunov vectors satisfying (2.18), where $^{\ast }$ denotes the complex conjugate. Any two non-collinear vectors, e.g. $\unicode[STIX]{x1D6FC}_{j}=\unicode[STIX]{x1D6FC}_{j}^{\ast }=1/2$ or $\unicode[STIX]{x1D6FC}_{j}=-\unicode[STIX]{x1D6FC}_{j}^{\ast }=-\text{i}/2$ , can be taken to define the Lyapunov subspace

(3.8) $$\begin{eqnarray}\left[\begin{array}{@{}cc@{}}\cos (\unicode[STIX]{x1D714}_{j}t) & -\sin (\unicode[STIX]{x1D714}_{j}t)\\ \sin (\unicode[STIX]{x1D714}_{j}t) & \cos (\unicode[STIX]{x1D714}_{j}t)\end{array}\right]\left[\begin{array}{@{}c@{}}\text{Re}(\hat{\boldsymbol{q}_{j}})\\ \text{Im}(\hat{\boldsymbol{q}_{j}})\end{array}\right].\end{eqnarray}$$

On a fixed point, we showed that the plane spanned by the covariant Lyapunov vectors does not change in time because the plane spanned by $\text{Re}(\hat{\boldsymbol{q}_{j}})$ and $\text{Im}(\hat{\boldsymbol{q}_{j}})$ is constant. In other words, while the angles between covariant Lyapunov vectors in the different Lyapunov subspaces vary in time, the angles between different Lyapunov subspaces are constant. This is in contrast to the chaotic case, where the angles between Lyapunov subspaces vary in time. Using the subspaces instead of the covariant Lyapunov vectors is crucial for the analysis of chaotic thermoacoustic systems, as shown in § 5.3.

3.2 Floquet analysis of limit cycles: connection with covariant Lyapunov vectors

Similarly to § 3.1, in this section we show that if the attractor is a limit cycle, the Lyapunov exponents correspond to the real part of the Floquet exponents and that the covariant Lyapunov vectors correspond to the eigenvectors of the monodromy matrix. Consider the tangent problem (2.3). We assume that the solution is a limit cycle, i.e. the solution is $T$ -periodic, i.e. $\bar{\boldsymbol{q}}(t+T)=\bar{\boldsymbol{q}}(t)$ , hence, the Jacobian is $T$ -periodic, i.e. $\unicode[STIX]{x1D645}\equiv \text{d}\boldsymbol{F}/\text{d}\boldsymbol{q}|_{\bar{\boldsymbol{q}}(t)}$ . Let $\unicode[STIX]{x1D64C}(t)=[\unicode[STIX]{x1D64C}_{1}|\cdots |\unicode[STIX]{x1D64C}_{N}]$ be the fundamental matrix and $\unicode[STIX]{x1D63D}$ the monodromy matrix (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983), i.e.

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}^{\prime }(t)=\unicode[STIX]{x1D64C}(t)\boldsymbol{c}, & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64C}(t+T)=\unicode[STIX]{x1D64C}(t)\unicode[STIX]{x1D63D}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{c}$ is the initial condition, $\boldsymbol{q}^{\prime }(0)$ , in the basis $\{\unicode[STIX]{x1D64C}_{1}(0),\ldots ,\unicode[STIX]{x1D64C}_{N}(0)\}$ . Let $\boldsymbol{b}_{j}$ be the eigenvector of $\unicode[STIX]{x1D63D}$ corresponding to the Floquet multiplier $\unicode[STIX]{x1D70C}_{j}=\text{e}^{\unicode[STIX]{x1D708}_{j}T}$ , where $\unicode[STIX]{x1D708}_{j}$ is the $j$ th Floquet exponent, sorted in descending order according to its real part, i.e. $j=1$ denotes the Floquet exponent with largest growth rate. Although the Floquet multipliers, which are the eigenvalues of the linearised Poincaré map, $\unicode[STIX]{x1D64C}(t)$ , are not unique because of the periodicity of the complex exponential, the Floquet exponents are unique. Noting that $\boldsymbol{q}^{\prime }(t)=\boldsymbol{q}^{\prime }(t^{+}+mT)=\unicode[STIX]{x1D64C}(t^{+})\unicode[STIX]{x1D63D}^{m}\boldsymbol{c}$ , with $0\leqslant t^{+}<T$ and $m=0,1,2,\ldots ,$ and decomposing $\boldsymbol{c}$ in the eigenbasis $\{\boldsymbol{b}_{1},\ldots ,\boldsymbol{b}_{N}\}$ , where we abuse notation and re-use the symbol $\unicode[STIX]{x1D6FC}_{j}$ to represent the coordinates in the local basis, yields

(3.11) $$\begin{eqnarray}\displaystyle \boldsymbol{q}^{\prime }(t) & = & \displaystyle \boldsymbol{q}^{\prime }(t^{+}+mT)\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D64C}(t^{+})\unicode[STIX]{x1D63D}^{m}\boldsymbol{c}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D64C}(t^{+})\unicode[STIX]{x1D63D}^{m}\mathop{\sum }_{j=1}^{N}(\underbrace{\boldsymbol{c}\boldsymbol{\cdot }\boldsymbol{b}_{j}}_{\triangleq \unicode[STIX]{x1D6FC}_{j}})\boldsymbol{b}_{j}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D64C}(t^{+})\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\unicode[STIX]{x1D70C}_{j}^{m}\boldsymbol{b}_{j}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D64C}(t^{+})\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\text{e}^{\unicode[STIX]{x1D708}_{j}mT}\boldsymbol{b}_{j}.\end{eqnarray}$$

Using (2.7), we can calculate the Lyapunov exponent restricted to times $t=t^{+}+mT$

(3.12)

which shows that the $k$ th Lyapunov exponent, $\unicode[STIX]{x1D706}_{k}$ , is equal to the real part of the $k$ th Floquet exponent, $\text{Re}(\unicode[STIX]{x1D708}_{k})$ . The result is independent of $t^{+}$ and valid for any $t^{+}\in [0,T)$ .

Consider a covariant Lyapunov vector $\unicode[STIX]{x1D753}$ associated with the Lyapunov exponent $\unicode[STIX]{x1D706}=\text{Re}(\unicode[STIX]{x1D708})$ . A priori, we do not know the shape of $\unicode[STIX]{x1D753}$ . Notwithstanding, we can express it in terms of a slightly modified eigenbasis of $\unicode[STIX]{x1D63D}$

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D753}(t)=\unicode[STIX]{x1D64C}(t)\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}\text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})t}.\end{eqnarray}$$

By differentiating (3.13) in time, we obtain

(3.14) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D753}} & = & \displaystyle \dot{\unicode[STIX]{x1D64C}}\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}\text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})t}-\unicode[STIX]{x1D64C}\mathop{\sum }_{j=1}^{N}\text{Re}(\unicode[STIX]{x1D708}_{j})\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}\text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})t}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\unicode[STIX]{x1D64C}\mathop{\sum }_{j=1}^{N}\text{Re}(\unicode[STIX]{x1D708}_{j})\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}\text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})t}.\end{eqnarray}$$

Similarly to § 3.1, we consider linear combinations of modes that have the same value of Lyapunov exponent, i.e. $\unicode[STIX]{x1D6FC}_{k}=0:\unicode[STIX]{x1D706}_{k}\neq \unicode[STIX]{x1D706}$ . Thus

(3.15) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D753}} & = & \displaystyle \unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\unicode[STIX]{x1D64C}\mathop{\sum }_{j=1}^{N}\text{Re}(\unicode[STIX]{x1D708}_{j})\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}\text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})t}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\text{Re}(\unicode[STIX]{x1D708})\unicode[STIX]{x1D64C}\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}\text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})t}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D645}\unicode[STIX]{x1D753}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D753},\end{eqnarray}$$

which recovers the covariant Lyapunov vector equation (2.18). In the same vein as in the fixed-point case (§ 3.1), any real linear combination of Floquet modes that have the same growth rate is a covariant Lyapunov vector. Finally, notice that $\unicode[STIX]{x1D753}(t)$ need not be periodic, except if it spans a one-dimensional Lyapunov subspace. In fact, if it is a linear combination of two complex conjugate Floquet modes, we have

(3.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D753}(t+T) & = & \displaystyle \text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})(t+T)}\unicode[STIX]{x1D64C}(t+T)(\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}+\unicode[STIX]{x1D6FC}_{j}^{\ast }\boldsymbol{b}_{j}^{\ast })\nonumber\\ \displaystyle & = & \displaystyle \text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})(t+T)}\unicode[STIX]{x1D64C}(t)\unicode[STIX]{x1D63D}(\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}+\unicode[STIX]{x1D6FC}_{j}^{\ast }\boldsymbol{b}_{j}^{\ast })\nonumber\\ \displaystyle & = & \displaystyle \text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})(t+T)}\unicode[STIX]{x1D64C}(t)(\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}\text{e}^{\unicode[STIX]{x1D708}_{j}T}+\unicode[STIX]{x1D6FC}_{j}^{\ast }\boldsymbol{b}_{j}^{\ast }\text{e}^{\unicode[STIX]{x1D708}_{j}^{\ast }T})\nonumber\\ \displaystyle & = & \displaystyle \text{e}^{-\text{Re}(\unicode[STIX]{x1D708}_{j})t}\unicode[STIX]{x1D64C}(t)(\unicode[STIX]{x1D6FC}_{j}\boldsymbol{b}_{j}\text{e}^{\text{Im}(\unicode[STIX]{x1D708}_{j})T}+\unicode[STIX]{x1D6FC}_{j}^{\ast }\boldsymbol{b}_{j}^{\ast }\text{e}^{\text{Im}(\unicode[STIX]{x1D708}_{j}^{\ast })T})\nonumber\\ \displaystyle & \neq & \displaystyle \unicode[STIX]{x1D753}(t),\end{eqnarray}$$

which shows that $\unicode[STIX]{x1D753}$ is not $T$ -periodic. Although this result might seem odd at first, $\unicode[STIX]{x1D753}$ behaves similarly to the covariant Lyapunov vectors in (3.8) because, in both cases, the imaginary part dictates the rate at which they rotate in the plane spanned by the corresponding mode. Although the mathematics is more involved, the connection between Floquet analysis and covariant Lyapunov vector analysis naturally follows the connection with eigenvalue analysis of fixed points (§ 3.1): a limit cycle can be viewed as a fixed point of a Poincaré map. In summary, on the one hand, covariant Lyapunov vector analysis provides the same linear dynamics as eigenvalue (Floquet) analysis when $\bar{\boldsymbol{q}}$ is a fixed-point (periodic solution). On the other hand, covariant Lyapunov vector analysis provides the linear dynamics when $\bar{\boldsymbol{q}}$ is a chaotic attractor, where eigenvalue and Floquet analyses cannot be applied.

The general theoretical analysis we have presented, which can be applied to other problems in flow stability, concludes the first part of this paper. From now on, we focus on a chaotic thermoacoustic system, which is a multi-physical problem in thermo-fluid dynamics that is relevant to aeronautical propulsion and clean power generation.

4.1 Acoustics and heat source

A thermoacoustic system consists of three subsystems that interact with each other: the acoustics, flame and hydrodynamics (Lieuwen Reference Lieuwen2012; Magri Reference Magri2018). The acoustics strongly depend on the geometry of the configuration and the boundary conditions. The flame is governed by chemistry mechanisms and their interaction with the turbulent environment. The hydrodynamics is governed by the geometry of the inlets and flame holders, which generate large coherent structures due to flow instabilities (vortex shedding, shear layer instabilities, etc.), which, in turn, are modulated by turbulence. To accurately model thermoacoustic instabilities, high-fidelity simulations can be employed (e.g. Poinsot Reference Poinsot2017). However, in this fundamental paper, we aim at capturing the essential physical mechanisms of chaotic thermoacoustic instabilities. Therefore, we choose a prototypical time-delayed thermoacoustic system with a longitudinal acoustic cavity and a heat source modelled with a nonlinear time-delayed model (Subramanian et al. Reference Subramanian, Mariappan, Sujith and Wahi2011). The main assumptions are: (i) the acoustics are small perturbations onto a mean flow at rest with uniform density; (ii) viscosity and diffusivity are negligible; and (iii) the acoustics are one-dimensional, i.e. the cut-on frequency of the duct is much higher than the frequency of the instability. Under these assumptions, the linearisation of the inviscid momentum and energy equations yields, respectively (Balasubramanian & Sujith Reference Balasubramanian and Sujith2008; Juniper Reference Juniper2011; Magri & Juniper Reference Magri and Juniper2013)

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}=0, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D701}p-\dot{q}\unicode[STIX]{x1D6FF}(x-x_{f})=0, & \displaystyle\end{eqnarray}$$

where $u$ , $p$ , $\dot{q}$ , $x$ and $t$ are the non-dimensional velocity, pressure, heat-release rate, axial coordinate and time, respectively. The reference scales for speed, pressure, length and time are the mean-flow convection velocity, the mean-flow Mach number multiplied by the heat capacity factor, the length of the tube and the length of the tube divided by the mean-flow speed of sound, respectively. The duct has open ends, which means that the acoustic pressure is zero at the boundaries. The damping coefficient, $\unicode[STIX]{x1D701}$ , takes into account all the acoustic dissipation (§ 4.1.1). The spatial extent of the heat source is assumed to negligible as compared to the acoustic wavelength (Dowling Reference Dowling1997), thus, it is modelled as a compact source of acoustic energy through a Dirac delta (generalised) function, $\unicode[STIX]{x1D6FF}(x-x_{f})$ localised at $x_{f}=0.2$ . The heat-release rate is provided by a modified King’s law (King Reference King1914; Heckl Reference Heckl1988, Reference Heckl1990; Polifke et al. Reference Polifke, Poncet, Paschereit and Döbbeling2001; Orchini, Rigas & Juniper Reference Orchini, Rigas and Juniper2016)

(4.3) $$\begin{eqnarray}\dot{q}(t)=\unicode[STIX]{x1D6FD}[(1+u_{f}(t-\unicode[STIX]{x1D70F}))^{1/2}-1],\end{eqnarray}$$

which is a nonlinear time-delayed model for an electrically heated mesh of wires. This model has similar features to flame models, such as the $n$ – $\unicode[STIX]{x1D70F}$ model (e.g. Juniper & Sujith Reference Juniper and Sujith2018). In future work, the heat-release rate, $\dot{q}(t)$ , can be obtained, for example, from the dynamics of premixed flames (Kashinath, Hemchandra & Juniper Reference Kashinath, Hemchandra and Juniper2013a ,Reference Kashinath, Hemchandra and Juniper b ; Waugh et al. Reference Waugh, Kashinath and Juniper2014; Orchini et al. Reference Orchini, Illingworth and Juniper2015; Yu, Juniper & Magri Reference Yu, Juniper and Magri2019) or diffusion flames (Tyagi, Jamadar & Chakravarthy Reference Tyagi, Jamadar and Chakravarthy2007; Magri & Juniper Reference Magri and Juniper2014). Solving for the flame dynamics adds many numerical degrees of freedom to the state vector, which makes the calculations computationally more expensive, but it does not change the framework we propose. Because Lyapunov analysis is valid only for smooth dynamical systems, we approximate the heat-release law (4.3) by a fourth-degree polynomial in a small neighbourhood of $u_{f}(t-\unicode[STIX]{x1D70F})=-1$ to make the derivative smooth (figure 2).

Figure 2. Comparison between King’s law and our smoothed version. The smoothed version is exactly equal to King’s law outside the range $|1+u_{f}|<\unicode[STIX]{x1D716}$ , inside of which it is approximated by a fourth-degree polynomial, which enables continuity of both the function and its derivative.

The heat parameter, $\unicode[STIX]{x1D6FD}$ , and time delay, $\unicode[STIX]{x1D70F}$ , encapsulate all information about the heat source, base velocity and ambient conditions. We transform the time-delayed problem into an initial value problem. (This operation is not mandatory, however, it enables us to use the adaptive initial value problem time integrator scipy.integrate.odeint with no modification.) Thus, we model the advection of a dummy variable $v$ with velocity $\unicode[STIX]{x1D70F}^{-1}$ as

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D70F}}\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}X}=0,\quad 0\leqslant X\leqslant 1, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle v(X=0,t)=u_{f}(t). & \displaystyle\end{eqnarray}$$

The dummy variable $v$ takes time $\unicode[STIX]{x1D70F}$ to travel from the left to the right boundary. Therefore, the time-delayed acoustic velocity is provided by the value of $v$ at the right boundary, i.e. $u_{f}(t-\unicode[STIX]{x1D70F})=v(X=1,t)$ . The calculation of the time-delayed acoustic velocity via (4.4) adds only a few numerical degrees of freedom (§ 4.1.1). The time-delayed problem (4.1)–(4.3) is mathematically equivalent to the initial value problem (4.1)–(4.2) and (4.4)–(4.5) (Jarlebring Reference Jarlebring2008).

4.1.1 Numerical discretisation

Equations (4.1), (4.2) are discretised by a Galerkin method (Zinn & Lores Reference Zinn and Lores1971). First, the acoustic variables are separated in time and space as

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle u(x,t)=\mathop{\sum }_{j=1}^{N_{g}}\unicode[STIX]{x1D702}_{j}(t)\cos (j\unicode[STIX]{x03C0}x), & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle p(x,t)=-\mathop{\sum }_{j=1}^{N_{g}}\unicode[STIX]{x1D707}_{j}(t)\sin (j\unicode[STIX]{x03C0}x), & \displaystyle\end{eqnarray}$$

where each spatial function is a natural acoustic mode of the open-ended duct. The partial differential equations (4.1), (4.2) are projected onto the Galerkin spatial basis $\{\cos (\unicode[STIX]{x03C0}x),\cos (2\unicode[STIX]{x03C0}x),\ldots ,\cos (N_{g}\unicode[STIX]{x03C0}x)\}$ to yield

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D702}}_{j}-j\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{j}=0, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D707}}_{j}+j\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{j}+\unicode[STIX]{x1D701}_{j}\unicode[STIX]{x1D707}_{j}+2\dot{q}\sin (j\unicode[STIX]{x03C0}x_{f})=0. & \displaystyle\end{eqnarray}$$

The system has $2N_{g}$ degrees of freedom. The time-delayed velocity becomes

(4.10) $$\begin{eqnarray}u_{f}(t-\unicode[STIX]{x1D70F})=\mathop{\sum }_{k=1}^{N_{g}}\unicode[STIX]{x1D702}_{k}(t-\unicode[STIX]{x1D70F})\cos (k\unicode[STIX]{x03C0}x_{f}),\end{eqnarray}$$

and the damping, $\unicode[STIX]{x1D701}_{j}$ , is modelled by a modal expression that damps out higher-frequency oscillations, $\unicode[STIX]{x1D701}_{j}=c_{1}\,j^{2}+c_{2}\,j^{1/2}$ , where $c_{1}=0.1$ and $c_{2}=0.06$ (Subramanian et al. Reference Subramanian, Mariappan, Sujith and Wahi2011). This damping model originates from physical principles, as explained in Landau & Lifshitz (Reference Landau and Lifshitz1987). We have assumed that the mean flow is sufficiently slow such that it can be neglected. Adding a mean flow may quantitatively change the phases between acoustic waves (Dowling & Morgans Reference Dowling and Morgans2005), but the conclusions of this paper are qualitatively unaffected. With a mean flow, a wave approach can be used instead (Dowling & Morgans Reference Dowling and Morgans2005).

The linear advection equation (4.4) is discretised using $N_{c}+1$ points with a Chebyshev spectral method (Trefethen Reference Trefethen2000). This discretisation adds $N_{c}$ degrees of freedom. The resulting discretised system is time integrated using scipy.integrate. odeint, which in turn calls lsoda from the odepack library. This method detects numerical stiffness and switches automatically between the Adams method in the non-stiff case and a backward differentiation formula in the stiff case (Petzold Reference Petzold1983). On numerical discretisation, the thermoacoustic state vector is $\bar{\boldsymbol{q}}$ , with $\bar{\boldsymbol{q}}^{\text{T}}=$ $(\unicode[STIX]{x1D702}_{1},\ldots ,\unicode[STIX]{x1D702}_{N_{g}},\unicode[STIX]{x1D707}_{1},\ldots ,\unicode[STIX]{x1D707}_{N_{g}},v_{1},\ldots ,v_{N_{c}})$ .

4.1.2 Effect of numerical discretisation

To investigate the effect of the numerical discretisation, we perform two types of tests. Figure 3 depicts the first type of test, which is done on a chaotic attractor. Figure 3(a) shows the standard deviation of $\unicode[STIX]{x1D702}_{j}(t),j=1,\ldots ,N_{g}$ for different values of $N_{g}$ with fixed $N_{c}=10$ . The dominant unstable mode is the first mode because the heat source is located at $x_{f}=0.2$ , where most of the energy excites the first mode. Apart from $N_{g}=5$ , where a large difference is observed, the calculations with small $N_{g}$ correctly capture the energy associated with each of the Galerkin modes that they compute. When increasing the number of Galerkin modes, the accuracy on the modes that were previously included does not improve. The benefit of increasing $N_{g}$ is to increase the spatial resolution by including higher wavenumbers. The magnitude of the standard deviation decays sharply at the beginning up to $j=10$ , followed by a slower decay. Therefore, capturing modes of lower intensity, e.g. $O(\text{SD}[\unicode[STIX]{x1D702}_{j}])\sim 10^{-3}$ , requires a large increase in the number of Galerkin modes. Thus, a good compromise between accuracy and computational cost is obtained with $N_{g}=10$ , which is the number of Galerkin modes used throughout the rest of this paper.

Figure 3. Convergence study on a chaotic attractor ( $\unicode[STIX]{x1D6FD}=7.0$ , $\unicode[STIX]{x1D70F}=0.2$ ). (a) Standard deviation of the Galerkin modes $\unicode[STIX]{x1D702}_{j}(t)$ on a chaotic attractor, varying the number of Galerkin modes, $N_{g}$ , with the number of Chebyshev points fixed to $N_{c}+1=11$ . The fifth and tenth modes are of the order of machine precision. (b) Same as panel (a), but with the number of Galerkin modes fixed to $N_{g}=10$ and number of Chebyshev points, $N_{c}+1$ , varying.

In § 3.1, we showed analytically that, if the attractor is a fixed point, the Lyapunov exponents are equal to the real part of the eigenvalues of the Jacobian at the fixed point. Therefore, the difference between the two is a metric that can be used to assess the quality of the numerical solution, which is the second test. Figure 4(a) shows the evolution of the first four Lyapunov exponents, the real part of the corresponding eigenvalues and their converged values as the number of Galerkin modes is increased to 30. As explained in § 5, the Lyapunov spectrum of this system is composed of Lyapunov exponents with double multiplicity. The residuals of these Lyapunov exponents are of the order of $10^{-4}$ at the end of the simulations. Once again, a good compromise between the accuracy of the first four Lyapunov exponents and computational cost is obtained with $N_{g}=10$ . A similar analysis is run by fixing the Galerkin modes to $N_{g}=10$ and varying the Chebyshev points, $N_{c}+1$ , (figures 3 b, 4 b), which shows that the influence of the number of Chebyshev points is not significant. In this paper, we use $N_{c}+1=11$ .

Figure 4. Convergence study on a fixed point ( $\unicode[STIX]{x1D6FD}=0.4$ , $\unicode[STIX]{x1D70F}=0.2$ ). (a) First four Lyapunov exponents (since they come in pairs, their mean is plotted) and real part of the corresponding eigenvalues, varying the number of Galerkin modes, $N_{g}$ , with the number of Chebyshev points fixed to $N_{c}+1=11$ . The vertical axes’ ranges correspond to $\pm 0.5\,\%$ of the converged value (dotted line). (b) Same as panel (a), but with the number of Galerkin modes fixed to $N_{g}=10$ and number of Chebyshev points, $N_{c}+1$ , varying. The vertical axes’ ranges correspond to $\pm 5\,\%$ of the converged value (dotted line).

5.1 Analysis on a fixed-point attractor

By setting $\unicode[STIX]{x1D6FD}=0.4$ , the solution converges to the fixed point $\bar{\boldsymbol{q}}=0$ . The spinup and spindown times are the same and equal to 200 time units, while the simulation time is 1000 with a time segment of 0.01, yielding $[200,800]$ as the interval of study. The spinup time is chosen a posteriori, by choosing a time such that the system is past the transient regime. Figure 7 shows the Lyapunov exponents and the real part of the eigenvalues of the Jacobian of the system at $\bar{\boldsymbol{q}}=0$ , demonstrating that the Lyapunov spectrum matches the real part of the eigenvalues. There are 16 distinct values of Lyapunov exponents, 14 of which have multiplicity two, corresponding to Lyapunov subspaces of dimension 2, while $\unicode[STIX]{x1D706}_{25}$ and $\unicode[STIX]{x1D706}_{30}$ correspond to one-dimensional Lyapunov subspaces, as described in § 3.1, for a total of 16 Lyapunov subspaces: $\unicode[STIX]{x1D6FA}_{1},\ldots ,\unicode[STIX]{x1D6FA}_{16}$ . The velocity and pressure components of the first Galerkin mode, $\unicode[STIX]{x1D702}_{1}$ and $\unicode[STIX]{x1D707}_{1}$ , of the covariant Lyapunov vectors $\unicode[STIX]{x1D753}^{(1)}$ and $\unicode[STIX]{x1D753}^{(2)}$ are denoted $\unicode[STIX]{x1D753}_{1}^{(1)}$ , $\unicode[STIX]{x1D753}_{11}^{(1)}$ , $\unicode[STIX]{x1D753}_{1}^{(2)}$ and $\unicode[STIX]{x1D753}_{11}^{(2)}$ , respectively (note that the 11th component of the state vector corresponds to $\unicode[STIX]{x1D707}_{1}$ in the arrangement of § 4.1.1). They are plotted in figure 8(a). The time series are not purely sinusoidal, as predicted by the analytical result (3.8) of § 3.1, because the covariant Lyapunov vectors are defined up to a time-varying factor. This time-varying factor is the normalisation that is imposed in the Gram–Schmidt orthonormalisation (i.e. the QR decomposition) in step (iv) in § 2.3. This normalisation varies in time because it is repeated at every time segment. Therefore, the time-varying normalisation generates higher harmonics in the power spectral density. This can be seen by comparing the power spectral densities of $\unicode[STIX]{x1D753}_{1}^{(1)}$ , $V_{1}$ and $V_{1}/\Vert \boldsymbol{V}\Vert$ (figure 8 b), where $\boldsymbol{V}$ is the first vector of (3.8). While $V_{1}$ expectedly presents one mode only at $f_{1}=\unicode[STIX]{x1D714}_{1}/2\unicode[STIX]{x03C0}\approx 0.6$ , $V_{1}/\Vert \boldsymbol{V}\Vert$ has peaks at frequencies of the form $kf_{1},k\in \{1,3,5,\ldots \}$ , exactly like $\unicode[STIX]{x1D753}_{1}^{(1)}$ . The mean of the angles between the Lyapunov subspaces and the eigensubspaces are shown in figure 9. (Eigensubspaces are subspaces that are spanned by eigenvectors corresponding to eigenvalues with the same real part, e.g. a pair of complex conjugates, as described in § 3.1.) The main diagonal, which compares Lyapunov subspaces to eigensubspaces of the same growth rate, is 0 (to precision), showing that the Lyapunov subspaces are indeed equal to the (constant) eigensubspaces. The analysis of this section numerically shows the equivalence between eigenvectors and covariant Lyapunov vectors on stable fixed points, as analytically explained in § 3.1.

Figure 7. Real part of the eigenvalues (▪, blue) and Lyapunov spectrum (▪, navy blue) match on the fixed-point solution $\bar{\boldsymbol{q}}=0$ ( $\unicode[STIX]{x1D6FD}=0.4$ ).

Figure 8. Fixed-point solution $\bar{\boldsymbol{q}}=0$ ( $\unicode[STIX]{x1D6FD}=0.4$ ). (a) Velocity and pressure components of the first Galerkin mode, $\unicode[STIX]{x1D702}_{1}$ and $\unicode[STIX]{x1D707}_{1}$ , of the covariant Lyapunov vectors $\unicode[STIX]{x1D753}^{(1)}$ and $\unicode[STIX]{x1D753}^{(2)}$ : $\unicode[STIX]{x1D753}_{1}^{(1)}$ (——, blue), $\unicode[STIX]{x1D753}_{11}^{(1)}$ (——, red), $\unicode[STIX]{x1D753}_{1}^{(2)}$ (——, green) and $\unicode[STIX]{x1D753}_{11}^{(2)}$ (——, purple) versus time – each component oscillates between the corresponding component of $\pm \sqrt{\text{Re}(\hat{\boldsymbol{q}})^{2}+\text{Im}(\hat{\boldsymbol{q}})^{2}}$ , where $\hat{\boldsymbol{q}}$ is the corresponding eigenvector. (b) Power spectral density of $\unicode[STIX]{x1D753}_{1}^{(1)}$ (——, blue), $V_{1}/\Vert \boldsymbol{V}\Vert$ (——, red) and $V_{1}$ (——, green), where $\boldsymbol{V}$ is the first vector of (3.8). The vertical axis is logarithmic of base 10.

Figure 9. Mean angle between Lyapunov subspaces and eigensubspaces, $\angle (\unicode[STIX]{x1D6FA}_{a},{\mathcal{Q}}_{b})$ , on the fixed-point solution $\bar{\boldsymbol{q}}=0$ ( $\unicode[STIX]{x1D6FD}=0.4$ ). (Standard deviation not shown because it is 0 to precision.) The main diagonal, which corresponds to the angles between Lyapunov subspaces and eigenspaces of the same growth rate (same index), is 0, showing that covariant Lyapunov vectors or eigenvectors of the same growth rate span the same subspaces.

5.2 Analysis on a periodic attractor

By increasing the heat-release intensity parameter to $\unicode[STIX]{x1D6FD}=2.5$ , the thermoacoustic system converges to a limit cycle. The spinup and spindown times are the same and equal to 200 time units, while the simulation time is 1000 with a time segment of 0.01, yielding $[200,800]$ as the interval of study. The velocity at the heat source $u_{f}(t)$ (figure 10 a) oscillates within $[-4.13,4.73]$ . The fact that the minima and maxima of $u_{f}(t)$ are not equal in absolute value can be explained by the asymmetry of the heat-release law (4.3). The period is $T_{0}=1.95$ , corresponding to a frequency of $f_{0}=0.51$ , which appears in the power spectral density of $u_{f}(t)$ (figure 10 b) as a maximum peak. The subsequent peaks occur at $kf_{0},\,k\in \{2,3,\ldots \}$ . The frequency $f_{0}=0.51$ has a value that is close to the natural acoustic frequency of the first mode of the duct, which is $f=\unicode[STIX]{x03C0}/(2\unicode[STIX]{x03C0})=0.5$ (§ 4.1.1). The frequency shift is physically due to the effect of the heat release and damping.

Figure 10. Acoustic velocity at the heat-source location, $u_{f}(t)$ , of a limit cycle ( $\unicode[STIX]{x1D6FD}=2.5$ ). (a) Time series with initial transient (the full simulation period (1000) is not depicted). Here, $u_{f}(t)$ becomes periodic with period $1.95$ at $t\approx 30$ . (b) Power spectral density. The vertical axis is logarithmic of base 10. The global maximum is at $f_{0}=0.51$ and the other local maxima are its higher harmonics.

Figure 11. Real part of the Floquet exponents (▪, blue) and Lyapunov spectrum (▪, navy blue) (first 20) match on a limit cycle solution ( $\unicode[STIX]{x1D6FD}=2.5$ ). The remaining 10 Floquet exponents are not shown because their values are large in absolute value and the accuracy of their calculation is limited by machine precision in the computation of the monodromy matrix.

Figure 11 shows the real part of the first 20 Floquet exponent and the corresponding Lyapunov exponents. The remaining 10 Floquet exponents are not shown because their values are large in absolute value and the accuracy of their calculation is limited by machine precision in the computation of the monodromy matrix (e.g. $\text{e}^{\unicode[STIX]{x1D706}_{21}T}\approx \text{e}^{-20\times 2}\sim 10^{-18}$ ). These modes are non-physical and correspond to the Chebyshev discretisation of the advection equation, as described § 4.1.1. The first Lyapunov and Floquet exponents are zero (to a small numerical error), i.e. they are the neutral modes, and they correspond to vectors tangent to the limit cycle. The second Lyapunov exponent is the least stable mode, which corresponds to a one-dimensional Lyapunov subspace. The remaining Lyapunov exponents have multiplicity two and match the real part of the Floquet exponents, the latter of which come as pairs of complex conjugates. As explained in § 5.1 and the present section, the fact the Lyapunov subspaces have double multiplicity has a physical interpretation. The thermoacoustic dynamics is driven by the nonlinear saturation of the thermoacoustic eigenfunctions, which come as complex conjugate pairs. As shown in § 5.3, this physical mechanism is dominant even when the system is chaotic.

Figures 12(a), 12(b) show the mean and the standard deviation, respectively, of the angle between the Lyapunov subspaces and the Floquet subspaces (subspaces spanned by groups of eigenvectors of the monodromy matrix that have the same real part of the Floquet exponent). The fact that the main diagonal of figure 12(b) is 0 demonstrates that these angles are constant, while the fact that the main diagonal of figure 12(a) is also 0 shows that the constant angles are 0. Thus, the Lyapunov subspaces are equal to the Floquet subspaces. Similarly to § 5.1, the analysis of this section numerically shows the equivalence between Floquet vectors (eigenvectors of the monodromy matrix) and covariant Lyapunov vectors on stable limit cycles, as analytically explained in § 3.2.

Figure 12. Angles between Lyapunov subspaces and Floquet subspaces, $\angle (\unicode[STIX]{x1D6FA}_{a},{\mathcal{B}}_{b})$ , on a limit cycle solution ( $\unicode[STIX]{x1D6FD}=2.5$ ). (a) Mean. (b) Standard deviation. The main diagonals of both figures are 0, showing that Floquet vectors (eigenvectors of the monodromy matrix) or covariant Lyapunov vectors of the same growth rate span the same subspaces.

5.3 Analysis on a chaotic attractor

The heat-release intensity parameter is further increased to $\unicode[STIX]{x1D6FD}=7.0$ , with the system converging to a chaotic attractor. In a chaotic solution, only covariant Lyapunov vector analysis can calculate the linear dynamics of the attractor. Eigenvalue and Floquet analyses are no longer valid. The acoustic velocity at the base of the heat source, $u_{f}(t)$ , is oscillatory (figure 13 a) but aperiodic. While the dominant peak of the power spectral density (figure 13 b) is largely unaltered from § 5.2, the frequency spectrum is now denser and shows multiple peaks at several frequencies, which indicates the presence of chaos. However, quasi-periodic solutions can exhibit the same behaviour and can be mistaken for chaotic. It is hard thus to classify the attractor from figures 13(a), 13(b) alone. Instead, we use the Lyapunov spectrum (figure 14). The first Lyapunov exponent, $\unicode[STIX]{x1D706}_{1}\approx 0.12$ , is positive, thus confirming that the attractor is indeed chaotic; $\unicode[STIX]{x1D753}_{2}(t)$ is the neutral covariant Lyapunov vector because $\unicode[STIX]{x1D706}_{2}=0$ to numerical error. The remaining Lyapunov exponents correspond to one-dimensional modes, except for the pairs $(\unicode[STIX]{x1D706}_{9},\unicode[STIX]{x1D706}_{10})$ , $(\unicode[STIX]{x1D706}_{19},\unicode[STIX]{x1D706}_{20})$ , $(\unicode[STIX]{x1D706}_{21},\unicode[STIX]{x1D706}_{22})$ and $(\unicode[STIX]{x1D706}_{24},\unicode[STIX]{x1D706}_{25})$ , each corresponding to two-dimensional Lyapunov subspaces. Thus, we conclude that

(5.2) $$\begin{eqnarray}\displaystyle & \displaystyle E_{q}^{u}=\text{Span}(\unicode[STIX]{x1D753}_{1}), & \displaystyle\end{eqnarray}$$

(5.3) $$\begin{eqnarray}\displaystyle & \displaystyle E_{q}^{n}=\text{Span}(\unicode[STIX]{x1D753}_{2}), & \displaystyle\end{eqnarray}$$

(5.4) $$\begin{eqnarray}\displaystyle & \displaystyle E_{q}^{s}=\text{Span}(\unicode[STIX]{x1D753}_{3},\ldots ,\unicode[STIX]{x1D753}_{30}). & \displaystyle\end{eqnarray}$$

For sensitivities to exist, the angles between these subspaces must be bounded away from 0 (), that is, the attractor must be hyperbolic. In the next section, this question is investigated on multiple design points where the system exhibits chaotic behaviour.

Figure 13. Acoustic velocity at the heat-source location, $u_{f}(t)$ , on a limit cycle solution ( $\unicode[STIX]{x1D6FD}=2.5$ ). (a) Time series. $u_{f}(t)$ is oscillatory, but aperiodic. (b) Power spectral density. The global maximum is at $f_{0}=0.52$ , which is close to the global maximum found in the limit cycle case of § 5.2.

Figure 14. Complete Lyapunov spectrum, with a closeup of the first 3 Lyapunov exponents, on a chaotic attractor ( $\unicode[STIX]{x1D6FD}=7.0$ ).  $\unicode[STIX]{x1D706}_{1}\approx 0.12$ and $\unicode[STIX]{x1D706}_{2}=0$ to numerical error (neutral mode). The remaining Lyapunov exponents are negative.

5.4 Are chaotic acoustic attractors hyperbolic?

In § 6.2, the sensitivities of the time-averaged acoustic energy with respect to the heat-source parameters are embedded in an optimisation routine to minimise the size of acoustic oscillations. However, as discussed in § 2.6, for such sensitivities to exist, the thermoacoustic chaotic attractor should be hyperbolic (§ 2.4), otherwise the sensitivities of (1.1) might not exist. Here, we seek to determine the hyperbolicity of the system. To determine whether a system is hyperbolic, the complete spectrum should be computed to construct the unstable, neutral and stable subspaces. This is possible with the reduced-order model of this paper, but it could be prohibitively expensive in high-dimensional systems. For the latter, only a portion of the Lyapunov spectrum and covariant vectors is typically calculated (e.g. Blonigan et al. Reference Blonigan, Fernandez, Murman, Wang, Rigas and Magri2016; Fernandez & Wang Reference Fernandez and Wang2017). Because it is computationally expensive to determine hyperbolicity everywhere in the design space, we restrict ourselves to 9 points, roughly equally spaced inside the chaotic areas (green in figure 5), which are shown in figure 15(a) (labels A to I). The probability density functions (PDFs) of the angles between the three pairs of elements from $E^{u},E^{n},E^{s}$ are calculated at each of these design points (figure 15 b). On the one hand, design points D and E are not hyperbolic, since the PDF of $\unicode[STIX]{x1D703}_{n,s}$ is non-zero at $\unicode[STIX]{x1D703}=0$ , demonstrating that the system exhibits tangencies between these two subspaces. On the other hand, the PDFs of the remaining 7 points indicate that these are hyperbolic, which, physically means that time-averaged cost functionals respond smoothly to small changes in the design parameters, i.e. their sensitivities exist. (For completeness, we also report that there is evidence that shadowing-based methods work well in some non-hyperbolic systems (Ni Reference Ni2019).) In conclusion, we found that thermoacoustic systems can physically exhibit both hyperbolic and non-hyperbolic chaos, depending on the design point. Notwithstanding, starting from design point B, we will employ shadowing techniques to compute sensitivities, which are employed in an optimisation routine to minimise the acoustic energy (§ 6.2).

Figure 15. (a) Design points where the probability density functions are calculated. (b) Probability density function of angles between $E^{u}$ , $E^{n}$ , $E^{s}$ . Tangencies are observed for points D and E, demonstrating that the chaotic attractors in these positions are not hyperbolic.

6.1 Time-averaged cost functionals

We analyse the chaotic acoustic oscillations of point B shown in figure 15(a). Because thermoacoustics is a multi-physical phenomenon, there are different norms (Chu Reference Chu1965; George & Sujith Reference George and Sujith2012), semi-norms (Magri Reference Magri2015; Blumenthal et al. Reference Blumenthal, Tangirala, Sujith and Polifke2016) and functionals to define a physical measure. For thermoacoustic systems with negligible mean flow, which cannot advect flow inhomogeneities like entropy spots, the acoustic energy is a suitable quantity of interest. The instantaneous acoustic energy of the whole system is defined as

(6.1) $$\begin{eqnarray}\displaystyle E_{ac}(t)\triangleq {\textstyle \frac{1}{2}}\int _{0}^{1}(u^{2}(t)+p^{2}(t))\,\text{d}x, & & \displaystyle\end{eqnarray}$$

which is the sum of the acoustic kinetic and potential energies, i.e. it is the Hamiltonian (constant of motion) of the natural acoustic system. Because of Parseval’s theorem, the acoustic energy is related to the Galerkin modes as $E_{ac}(t)=\frac{1}{4}\sum _{j=1}^{N_{g}}(\unicode[STIX]{x1D702}_{j}^{2}(t)+\unicode[STIX]{x1D707}_{j}^{2}(t))$ . The acoustic energy, $E_{ac}$ , is (half) the Euclidean norm of the thermoacoustic system under investigation. In chaotic acoustic oscillations, we are interested in calculating the sensitivity of the time-averaged acoustic energy, $\langle E_{ac}\rangle$ . Figure 16 shows the acoustic energy in the refined area of figure 5, with the optimisation starting point, B, marked. Regions similar to those depicted by different colours in figure 15(a) are visible. Notably, a sharp discontinuity exists to the right of the chaotic regions, which corresponds to a line of bifurcation points. Furthermore, the time-averaged acoustic energy is multi-modal, exhibiting multiple local extrema. Interestingly, continuous regions of the same type of attractor (figure 6 b) do not have extrema in their interior. Instead, the extrema are found at the edges of the regions, which suggests that gradient-based optimisation algorithms will be capable of finding the boundaries that separate attractors.

Figure 16. Colour map of time-averaged acoustic energy, $\langle E_{ac}\rangle$ . To increase colour resolution near the starting point, B, points with $\langle E_{ac}\rangle \geqslant 30$ have the same colour.

6.2 Minimisation of acoustic energy in chaotic acoustic oscillations by optimal design

To minimise the acoustic energy, either the bifurcation diagram in the multidimensional parameter space is calculated (figure 16), which is computationally cumbersome, or a nonlinear optimisation problem of a time-averaged cost functional is solved, which is computationally affordable. Following the latter route, the optimisation problem is formulated as

(6.2) $$\begin{eqnarray}\left.\begin{array}{@{}ll@{}}\underset{\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D70F}}{\text{minimise}} & \left\langle E_{ac}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D70F})\right\rangle \\ \text{subject to} & \text{(4.1)},\text{(4.2)},\text{(4.3)},\text{(4.4)},\text{(4.5)},\text{(4.6)},\text{(4.7)}.\end{array}\right\}\end{eqnarray}$$

The set of parameters is updated via the sequential least squares programming method of the SciPy library. The optimisation stops when the condition

(6.3) $$\begin{eqnarray}\frac{\langle E_{ac}\rangle _{j}-\langle E_{ac}\rangle _{j+1}}{\langle E_{ac}\rangle _{j}}<0.01\end{eqnarray}$$

is met, that is, when the improvement between optimisation iterations is less than 1 % of the previous value. The usual gradient vanishing condition of extrema is not applied because of the existence of discontinuities. Since the gradient does not exist at these points, its numerical value cannot be trusted close to such points, which is why condition (6.3) is used instead. Figure 17(a) shows the cost functional as a function of the iteration in the optimisation algorithm. The acoustic energy decreases rapidly until iteration $7$ , where no progress is made and (6.3) is verified, indicating that the algorithm has converged to a local minimum ( $\unicode[STIX]{x1D6FD}=6.79889$ , $\unicode[STIX]{x1D70F}=0.18685$ ). Overall, the optimisation achieves a reduction of 14.8 % in acoustic energy in 7 iterations. Figure 17(b) shows the path that the optimisation procedure takes in the design space, showing that the final design is indeed a local minimum and that it is located at the edge of the region of chaotic solutions, as hypothesised in § 6.1. In conclusion, we found the set of parameters that produce a local minimum of the time-averaged acoustic energy of a chaotic thermoacoustic system. A similar algorithm can be used to find local maxima for maximal energy extraction in the design of thermoacoustic engines.

Figure 17. (a) Time-averaged acoustic energy versus optimisation iteration. The optimisation algorithm consecutively reduces the acoustic energy until iteration 7, where no progress is made. (b) Time-averaged acoustic energy, $\langle E_{ac}\rangle$ – colour map with optimisation path superimposed (——, red).

6.3 Future directions

In § 6.1, we found that the time-averaged acoustic energy, $\langle E_{ac}\rangle$ , displays an intricate behaviour: it is both discontinuous and multimodal. These two facts have strong implications when tackling an optimisation problem with a gradient-based approach. First, gradients are not defined at discontinuous points; second, gradient-based algorithms might converge to local minima, instead of the global optimum, as in § 6.2. While pure gradient-free algorithms might be computationally expensive, hybrid techniques, which couple gradient-based with gradient-free algorithms, can be a suitable compromise between overcoming the aforementioned issues of gradient-based algorithms and exploiting local gradient information to find the shortest path to an optimal design. For example, a Monte Carlo gradient-based optimisation, where multiple gradient-based optimisations are launched from various randomised initial design points and the best of the resulting (local) minima is chosen. Since thermoacoustic systems can exhibit rich dynamics and admit several types of solutions, a general technique for calculating sensitivities, including being capable of handling chaotic attractors, is required. As discussed in this work, covariant Lyapunov vector (CLV) analysis is the most general framework of analysis of these dynamical systems. Currently, shadowing methods, which exploit shadowing orbits via covariant Lyapunov vectors, are the leading candidate for sensitivity calculation in chaotic attractors.

For a practitioner, whose objective is to suppress oscillations, CLV techniques are required because chaos is always present in realistic turbulent systems. For example, in turbulent flows the hydrodynamic field modulates the heat released by the flame in a chaotic way (Lieuwen Reference Lieuwen2012). In the bistable region of a subcritical bifurcation, where the system is eigenvalue stable but a finite-amplitude periodic solution exists (Subramanian et al. Reference Subramanian, Mariappan, Sujith and Wahi2011), the turbulent hydrodynamic field chaotically forces the solution to oscillate around the stable fixed point or the upper branch limit cycle, depending on the initial condition. To eliminate the possibility of such chaotic limit cycles, the operating design point must be outside the hysteresis region. The boundary that separates these two regions (known as the fold point with one parameter) is marked by a discontinuity similar to the ones in figure 16. Because of the presence of chaos, it would not be possible to identify the fold point with traditional eigenvalue and Floquet analyses. However, it would be possible to identify it by using the CLV technique and optimisation we proposed (as shown in § 6.2).