A general framework for finding Hadamard instabilities

We return to the linear stability problem given in (3.5). A general procedure for detecting the presence or absence of Hadamard instabilities is developed. Since it is cast as an arbitrary matrix equation, there is no restriction on the dimensionality $N$ of the system, so our analysis in this subsection is applicable to models with any number of phases $n=N/2$ . Readers that would rather skip the linear algebra may proceed to the final paragraph of this subsection, where the method for determining posedness is recapitulated.

First, we bring (3.5) into a simpler form for analysis. The matrix $\unicode{x1D63C}$ must be invertible, in order for there to be $N$ independent time-evolving fields. Furthermore, we assume that the matrix $\unicode{x1D63C}^{-1}\unicode{x1D63F}$ is diagonalisable, since this covers all the specific cases in this paper. Then, the problem may be reformulated in terms of a basis $\{\hat {\boldsymbol{e}}_1,\ldots ,\hat {\boldsymbol{e}}_N\}$ with respect to which $\unicode{x1D63C}^{-1}\unicode{x1D63F}$ is diagonal. Therefore, for each matrix $\unicode{x1D648} \in \{\unicode{x1D63D}, \unicode{x1D63E}, \unicode{x1D63F}\}$ , we define

(4.1a,b) \begin{align} \kern2pt\hat{\kern-2pt \unicode{x1D648}} \,=\, \unicode{x1D64B}^{\,-1}\unicode{x1D63C}^{-1}\unicode{x1D648}\unicode{x1D64B}\quad \mathrm{and}\quad \hat {\boldsymbol{v}} = \unicode{x1D64B}^{\,-1}\boldsymbol{v}, \end{align}

for any vector $\boldsymbol{v}$ , where $\unicode{x1D64B}$ is a basis change matrix that diagonalises $\unicode{x1D63C}^{-1}\unicode{x1D63F}$ . With respect to this transformation, (3.5) becomes

(4.2) \begin{equation} \sigma \hat {\boldsymbol{r}} + \mathrm{i}\mkern 1mu k \kern2pt \hat{\kern-2pt \unicode{x1D63D}} \kern1pt\hat{\boldsymbol{r}} = \kern2pt \hat{\kern-2pt \unicode{x1D63E}} \kern1pt\hat {\boldsymbol{r}} - k^2 \kern2pt \hat{\kern-2pt \unicode{x1D63F}}\kern1pt \hat{\boldsymbol{r}}, \end{equation}

with $\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ a diagonal matrix. At high wavenumber $k\gg 1$ , we make the following asymptotic expansions:

(4.3a,b) \begin{equation} \sigma = -\sigma _2 k^2 - \mathrm{i}\mkern 1mu\sigma _1 k + \sigma _0 + \ldots , \quad \hat {\boldsymbol{r}} = \hat {\boldsymbol{r}}_0 + k^{-1}\hat {\boldsymbol{r}}_{-1} + \ldots , \end{equation}

substitute them into (4.2) and look for the leading-order terms. Therefore, at $O(k^2)$ , the problem reduces to

(4.4) \begin{equation} \kern2pt \hat{\kern-2pt \unicode{x1D63F}} \kern1pt \hat{\boldsymbol{r}}_0 = \sigma _2 \hat {\boldsymbol{r}}_0. \end{equation}

Noting the sign convention in (4.3a , b ), the eigenvalues $\sigma _2$ , which represent diffusion coefficients for the linear problem, must each have non-negative real part in order to avoid blow-up of $\textrm {Re}(\sigma )$ . The growth of modes with $\sigma _2 = 0$ is determined beyond this leading-order balance. If $\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ is not full rank, it has $i \in \{1,\ldots , N\}$ zero eigenvalues. Without loss of generality, we locate these in the first $i$ diagonal values of $\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ . The corresponding eigenvectors are determined only up to an $i$ -dimensional subspace ( $\hat {\boldsymbol{r}}_0 \in \textrm {span}\{\hat {\boldsymbol{e}}_1,\ldots ,\hat {\boldsymbol{e}}_i\}$ ), by (4.4).

Therefore, we proceed to the $O(k)$ part of the asymptotic expansion of (4.2). When $\sigma _2 = 0$ , this is

(4.5) \begin{equation} (\kern2pt \hat{\kern-2pt \unicode{x1D63D}} - \sigma _1 \unicode{x1D644}\,)\hat {\boldsymbol{r}}_0 = \mathrm{i}\mkern 1mu\kern2pt \hat{\kern-2pt \unicode{x1D63F}}\kern1pt \hat{\boldsymbol{r}}_{-1}. \end{equation}

Since $\hat {\boldsymbol{r}}_0 \in \textrm {span}\{\hat {\boldsymbol{e}}_1,\ldots ,\hat {\boldsymbol{e}}_i\}$ , only the first $i$ columns of $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}-\sigma _1\unicode{x1D644}$ enter into this system of equations on the left-hand side. Furthermore, only the first $i$ rows of (4.5) are needed to determine $\hat {\boldsymbol{r}}_0$ and these are rows for which the right-hand side is zero. Consequently, the $\sigma _1$ values are the eigenvalues of the matrix $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}$ with the last $N{-}i$ rows and columns removed. We shall write $\unicode{x1D648}_{\textit {red}}$ to denote any matrix $\unicode{x1D648}$ reduced in this way by deleting rows and columns associated with the nullspace of the diagonal matrix $\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ . Referring back to (4.3a , b ), we obtain a second criterion that must be met to avoid Hadamard instability: the eigenvalues $\sigma _1$ of $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ must be real. If these values are also distinct, then the growth rates stay bounded as $k \to \infty$ .

However, $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ may have repeated eigenvalues, which can also lead to blow-up of $\textrm {Re}(\sigma )$ . To see why, we proceed to the $O(1)$ equation with $\sigma _2 = 0$ , which reads

(4.6) \begin{equation} (\sigma _0 \unicode{x1D644} - \kern2pt \hat{\kern-2pt \unicode{x1D63E}}\,)\hat {\boldsymbol{r}}_0 + \mathrm{i}\mkern 1mu(\kern2pt \hat{\kern-2pt \unicode{x1D63D}} - \sigma _1 \unicode{x1D644}\,)\hat {\boldsymbol{r}}_{-1} = -\kern2pt \hat{\kern-2pt \unicode{x1D63F}}\kern1pt \hat{\boldsymbol{r}}_{-2}. \end{equation}

To eliminate dependence of the left-hand side on the unknown vector $\hat {\boldsymbol{r}}_{-2}$ , the left eigenvectors, corresponding to the eigenproblem adjoint to (4.2), may be used. By repeating the arguments used to determine $\hat {\boldsymbol{r}}_0$ , these may be expanded as $\hat {\boldsymbol{l}} = \hat {\boldsymbol{l}}_0 + k^{-1}\hat {\boldsymbol{l}}_{-1}+\ldots$ and inferred to satisfy $\hat {\boldsymbol{l}}_0^T\kern2pt \hat{\kern-2pt \unicode{x1D63F}}\,=\,\boldsymbol{0}$ and $\hat {\boldsymbol{l}}_0^T(\kern2pt \hat{\kern-2pt \unicode{x1D63D}} - \sigma _1 \unicode{x1D644}\,) = \mathrm{i}\mkern 1mu \hat {\boldsymbol{l}}_{-1}^T\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ (when $\sigma _2 = 0$ ). For any of the $i$ modes, the dot product of the leading-order left eigenvector $\hat {\boldsymbol{l}}_0$ may be taken with (4.6) and on rearranging the result, the formula

(4.7) \begin{equation} \sigma _0 = \frac {\hat {\boldsymbol{l}}_0 \boldsymbol{\cdot } \kern2pt \hat{\kern-2pt \unicode{x1D63E}}\,\hat {\boldsymbol{r}}_0 +\hat {\boldsymbol{l}}_{-1}\boldsymbol{\cdot }\kern2pt \hat{\kern-2pt \unicode{x1D63F}}\kern1pt \hat{\boldsymbol{r}}_{-1}}{\hat {\boldsymbol{l}}_0\boldsymbol{\cdot }\hat {\boldsymbol{r}}_0} \end{equation}

is obtained. Note that the relevant components of $\hat {\boldsymbol{l}}_{-1}$ and $\hat {\boldsymbol{r}}_{-1}$ required to compute the second term in the numerator are fully determined by inverting the final $N{-}i$ rows of (4.5) and their adjoint counterparts. The left and right eigenvectors for $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ are the vectors $\hat {\boldsymbol{l}}_0$ , $\hat {\boldsymbol{r}}_0$ with the last $N{-}i$ entries (which are all zeros) deleted. When $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ is diagonalisable, these vectors form a biorthonormal set, with the left and right eigenvectors for each mode satisfying $\hat {\boldsymbol{l}}_0\boldsymbol{\cdot }\hat {\boldsymbol{r}}_0 = 1$ , so the $O(1)$ growth rate in (4.7) is well defined. However, if $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ is not diagonalisable, at least one of its eigenvalues is defective. Therefore, $\sigma _1$ is a repeated eigenvalue associated with one or more Jordan chains of length at least two. Then for the full matrix $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}$ there are two pairs of corresponding generalised left and right eigenvectors $\hat {\boldsymbol{l}}_{0,1}$ , $\hat {\boldsymbol{l}}_{0,2}$ and $\hat {\boldsymbol{r}}_{0,1}$ , $\hat {\boldsymbol{r}}_{0,2}$ respectively (in $\textrm {span}\{\hat {\boldsymbol{e}}_1,\ldots ,\hat {\boldsymbol{e}}_i\}$ ), which satisfy

(4.8a,b) \begin{align} \begin{cases} \hat {\boldsymbol{l}}_{0,1}^T(\kern2pt \hat{\kern-2pt \unicode{x1D63D}}-\sigma _1 \unicode{x1D644}\,) = \mathrm{i}\mkern 1mu\hat {\boldsymbol{l}}_{-1}^T\kern2pt \hat{\kern-2pt \unicode{x1D63F}},\\[4pt] \hat {\boldsymbol{l}}_{0,2}^T(\kern2pt \hat{\kern-2pt \unicode{x1D63D}}-\sigma _1 \unicode{x1D644}\,) = \hat {\boldsymbol{l}}_{0,1}^T +\hat {\boldsymbol{\chi }}^T, \end{cases} \mathrm{and} \quad \begin{cases} (\kern2pt \hat{\kern-2pt \unicode{x1D63D}}-\sigma _1 \unicode{x1D644}\,)\hat {\boldsymbol{r}}_{0,1} = \mathrm{i}\mkern 1mu\kern2pt \hat{\kern-2pt \unicode{x1D63F}}\kern1pt \hat{\boldsymbol{r}}_{-1},\\[4pt] (\kern2pt \hat{\kern-2pt \unicode{x1D63D}}-\sigma _1 \unicode{x1D644}\,)\hat {\boldsymbol{r}}_{0,2} = \hat {\boldsymbol{r}}_{0,1} +\hat {\boldsymbol{\varGamma }}, \end{cases} \end{align}

where $\hat {\boldsymbol{r}}_{0,1} \equiv \hat {\boldsymbol{r}}_0$ and $\hat {\boldsymbol{l}}_{0,1}\equiv \hat {\boldsymbol{l}}_0$ , and $\hat {\boldsymbol{\chi }},\hat {\boldsymbol{\varGamma }}$ are unknown vectors in $\textrm{span} \{\hat{\boldsymbol{e}}_{i+1},\ldots ,\hat{\boldsymbol{e}}_N\}$ . In this case, the formula in (4.7) is always singular, since projecting any left eigenvector onto (4.8b ) shows that $\hat {\boldsymbol{l}}_0 \boldsymbol{\cdot }\hat {\boldsymbol{r}}_{0} = 0$ . Physically, this singularity can be thought to emerge from a resonance between two or more modes that collapse onto one another when $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ becomes defective. Examples of this are given below, in § 4.3.

The failure of (4.7) in these cases suggests the need for an alternative asymptotic expansion. Anticipating growth of some intermediate order between $O(k)$ and $O(1)$ , we replace the expansions in (4.3a , b ) with

(4.9a,b) \begin{equation} \sigma = -\mathrm{i}\mkern 1mu\sigma _1 k + \sigma _{1/2} k^{1/2} + \sigma _0 + \ldots , \quad \hat {\boldsymbol{r}} = \hat {\boldsymbol{r}}_{0,1} + k^{-1/2}\hat {\boldsymbol{r}}_{-1/2} + k^{-1}\hat {\boldsymbol{r}}_{-1}+\ldots . \end{equation}

This leaves the analysis at $O(k)$ unchanged and introduces the following equation at $O(k^{1/2})$ :

(4.10) \begin{equation} \sigma _{1/2} \hat {\boldsymbol{r}}_{0,1} + \mathrm{i}\mkern 1mu (\kern2pt \hat{\kern-2pt \unicode{x1D63D}} - \sigma _1 \unicode{x1D644}\,) \hat {\boldsymbol{r}}_{-1/2} + \kern2pt \hat{\kern-2pt \unicode{x1D63F}}\kern1pt \hat {\boldsymbol{r}}_{-3/2} = \boldsymbol{0}. \end{equation}

We project this onto $\hat {\boldsymbol{l}}_{0,2}$ and use (4.8), along with the fact that $\hat {\boldsymbol{l}}_{0,2}$ is orthogonal to the range of $\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ , to conclude that

(4.11) \begin{equation} \sigma _{1/2} \hat {\boldsymbol{l}}_{0,2}\boldsymbol{\cdot }\hat {\boldsymbol{r}}_{0,1} + \mathrm{i}\mkern 1mu(\hat {\boldsymbol{l}}_{0,1} +\hat {\boldsymbol{\chi }})\boldsymbol{\cdot }\hat {\boldsymbol{r}}_{-1/2} = 0. \end{equation}

The unknown vector $\hat {\boldsymbol{r}}_{-1/2}$ is eliminated by proceeding to the $O(1)$ equation. With the new expansion, this is

(4.12) \begin{equation} \sigma _{1/2}\hat {\boldsymbol{r}}_{-1/2} + \sigma _0 \hat {\boldsymbol{r}}_{0,1} +\mathrm{i}\mkern 1mu(\kern2pt \hat{\kern-2pt \unicode{x1D63D}} - \sigma _1 \unicode{x1D644}\,)\hat {\boldsymbol{r}}_{-1} - \hat {\unicode{x1D63E}}\,\hat {\boldsymbol{r}}_{0,1} + \kern2pt \hat{\kern-2pt \unicode{x1D63F}}\kern1pt \hat{\boldsymbol{r}}_{-2} = \boldsymbol{0}. \end{equation}

Then, we project this onto $\hat {\boldsymbol{l}}_{0,1}$ . Since $\hat {\boldsymbol{l}}_{0,1}\boldsymbol{\cdot }\hat {\boldsymbol{r}}_{0,1} = 0$ , the term containing $\sigma _0$ vanishes, along with the diffusive term which lies in an orthogonal subspace. After rearranging and using (4.11), as well as the $O(k^{3/2})$ part of the system, which implies that $\hat {\boldsymbol{r}}_{-1/2}\in \textrm {span}\{\hat {\boldsymbol{e}}_1,\ldots ,\hat {\boldsymbol{e}}_i\}$ , we obtain a formula for the $O(k^{1/2})$ part of the growth rate:

(4.13) \begin{equation} \sigma _{1/2} = \pm \frac {1 - \mathrm{i}\mkern 1mu}{2} \left ( \frac {2( \hat {\boldsymbol{l}}_{0,1}\boldsymbol{\cdot }\hat {\unicode{x1D63E}}\,\hat {\boldsymbol{r}}_{0,1} +\hat {\boldsymbol{l}}_{-1}\boldsymbol{\cdot }\kern2pt \hat{\kern-2pt \unicode{x1D63F}}\kern1pt \hat{\boldsymbol{r}}_{-1})} {\hat {\boldsymbol{l}}_{0,2} \boldsymbol{\cdot }\hat {\boldsymbol{r}}_{0,1} } \right )^{\!1/2}\!\!\!. \end{equation}

For Jordan chains of length two $\hat {\boldsymbol{l}}_{0,2}\boldsymbol{\cdot }\hat {\boldsymbol{r}}_{0,1}=|\hat {\boldsymbol{l}}_{0,2}||\hat {\boldsymbol{r}}_{0,1}|\neq 0$ , provided both the left and right vectors correspond to the same Jordan block. Consequently, (4.13) implies that there is a mode such that $\textrm {Re}(\sigma ) \sim k^{1/2}$ , provided the terms in the numerator do not interact in a way that causes it to vanish.

Conversely, for longer Jordan chains, the denominator in the (4.13) formula is also guaranteed to be singular. Different asymptotic expansions are needed, depending on the length of the the chain. However, to avoid these further complications, we terminate our analysis here, since cases where three or more modes intersect at high wavenumber are far less commonly encountered.

To summarise the analysis above, models up to second order that may be cast in the general form of (3.4) are ill posed as initial-value problems if any of the following conditions are met:

1. (i) Any eigenvalue of $\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ is negative, where $\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ denotes a diagonalisation of $\unicode{x1D63C}^{-1}\unicode{x1D63F}$ .

2. (ii) Any eigenvalue of $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ is complex, where $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ denotes the matrix formed by representing $\unicode{x1D63C}^{-1}\unicode{x1D63D}$ in the basis used to diagonalise $\unicode{x1D63C}^{-1}\unicode{x1D63F}$ in (i) and deleting each row and column $j$ such that the $j$ -th diagonal entry of $\kern2pt \hat{\kern-2pt \unicode{x1D63F}}$ is nonzero. We refer to $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ as a ‘reduced Jacobian’ in later analysis.

3. (iii) Repeated real eigenvalues of $\kern2pt \hat{\kern-2pt \unicode{x1D63D}}_{\textit {red}}$ of algebraic multiplicity $2$ share the same left and right eigenvectors $\hat {\boldsymbol{l}}_{0,1}$ and $\hat {\boldsymbol{r}}_{0,1}$ (up to normalisation), and the numerator of (4.13) is nonvanishing. (More generally, the expectation following from (4.7), is that repeated real eigenvalues of any algebraic multiplicity $m\geqslant 2$ imply ill posedness if the dimension of their associated eigenspace is strictly less than $m$ , but this is not explicitly proven above.)

For the remainder of this section, we apply these steps to different example systems.