Canonical form of structure constants and mutually orthogonal triad on the orbits of Bianchi cosmologies In Class. Quant. Grav. 28 (2011) 185007: "Linearization of homogeneous, nearly-isotropic cosmological models" at the start of section 2.3 the authors (Andrew Pontzen and Anthony Challinor) claim that we can simultaneously diagonalize the metric tensor restricted to the orbit, and bring the structure constants $C^k{}_{ij}$ into their canonical form. How can this be true in the general case?
Here we have
$$
[\xi_i,\xi_j] = -C^k{}_{ij}\xi_k,\tag{1}
$$
for translational killing vectors $\xi_i$ of some simply transitive 3-parameter symmetry group. The author previously states that we can decompose the structure constants as (this follows Landau-Lifshitz)
$$
C^k{}_{ij} = 2\delta^k_{[i}a_{j]} + \epsilon_{ij\ell}n^{\ell k},\tag{2}
$$
where $n^{ij} = n^{(ij)}$, and we achieve the canonical form by making a linear transformation of the form $\xi_i \mapsto \gamma_i^j\xi_j$ (letting $a,b,c,\ldots$ denote transformed indices) such that $n = \mathrm{diag}(n_1,n_2,n_3)$ and $a = (a, 0, 0)$. Finally one then normalizes $n$ and if possible $a$. By Landau-Lifshitz (page 112 in v.2) we first diagonalize $n^{ij}$ and then the Jacobi identity gives $n^{ij}a_j = 0$ so that $a_j$ is an eigenvector of $n^{ij}$ with eigenvalue 0, and we are free to set $a_j$ as above since we have imposed no further restriction on directions. Note that as far as I know, these are not tensorial indices; $a_i$ and $n^{ij}$ do not transform under the standard transformation laws.
Now, we have constructed some invariant frame field of the orbit (triad), defined by
$$
[e_i,\xi_j] = 0. \tag{3}
$$
The procedure outlined then first states that we can choose $e_i$ to be orthogonal (by selecting an orthogonal basis at some point, $p$, and Lie dragging it, by homogeneity of the metric with respect to $\xi_i$ and (3) it will be everywhere orthogonal on the orbit, my comment). Next we can choose $e_i|_p = \xi_i|_p$ (since the killing vectors span the tangent space at each point, my comment). Then follows the claim (note that we have $[e_i,e_j] = +C^k{}_{ij}e_k$ by this construction).

By making further linear reparametrizations of the $e_i$ and $\xi_i$, the $C^k{}_{ij}$ are brought into canonical form without disturbing the orthogonality.

This requires a transformation of the type $\xi_i \mapsto \gamma^j_i\xi_j$ that diagonalizes both the induced inner product $q_{ij}$ of $T_pH = T_pM|_H$ (where $H$ denotes the orbit) and $n^{ij}$. However, this is in general not possible for linear operators. Furthermore, once mutual diagonalization is achieved we are free to set $a = \mathrm{diag}(a,0,0)$ as in the canonical case if and only if for any two eigenvectors of $n^{ij}$ with eigenvalue $0$ the corresponding eigenvalues of $q_{ij}$ coincide.
 A: First off, we need to find how the matrix $n^{ij}$ transforms under the linear transformations considered. To this purpose note that we can construct the decomposition (2) in the OP by considering (following the convention the authors used)
$$
C^{ij} = C^i{}_{k\ell}\epsilon^{jk\ell}, \tag{4}
$$
which obeys $\epsilon_{jmn}C^{ij} = C^i{}_{k\ell}\epsilon^{k\ell}_{mn} = C^i{}_{[mn]} = C^i{}_{mn}$. Then we split $C^{ij} = C^{(ij)} + C^{[ij]}$ and take
\begin{align}\begin{split}
C^{(ij)} &\equiv n^{ij}, \\
C^{[ij]} &\equiv -\epsilon^{ijk}a_k,
\end{split}\tag{5}\end{align}
to define $n^{ij}$ and $a_k$. Here we used the fact that it is a three-parameter algebra, so that the anti-symmetric part and the vector have the same number of independent components. In other words
\begin{align}\begin{split}
n^{ij} &= C^{(i}{}_{k\ell}\epsilon^{j)k\ell}, \\
a_k &= -C^i{}_{ki}.
\end{split}\tag{5.1}\end{align}
Proceeding we will let $a,b,c$ denote transformed indices: $\xi_a = \gamma_a^j\xi_j$. We will also use $\gamma_i^a$ to denote the inverse of $\gamma^i_a$, letting the indices inform us which transformation we are dealing with. From (5.1) it is obvious that $a_c = \gamma_c^ja_j$, but the transformation of $n^{ij}$ is slightly trickier:
\begin{align}
n^{cd} &= C^k{}_{ij}\gamma^i_a\gamma^j_b\gamma^{(c}_k\epsilon^{d)ab} \\
&= \left(2\delta^k_{[i}a_{j]} + \epsilon_{ij\ell}n^{\ell k}\right)\gamma^i_a\gamma^j_b\gamma^{(c}_k\epsilon^{d)ab}. \tag{6}
\end{align}
Straight-forward calculation verifies that the first term in (6) vanishes and we have
\begin{align}
n^{cd} &= \epsilon^{ab(d}_{ij\ell}\gamma^{c)}_k\gamma^i_a\gamma^j_bn^{\ell k} \\
&= \epsilon^{abe}_{ij\ell}\gamma^{(d}_m\gamma^{c)}_k\gamma^m_e\gamma^i_a\gamma^j_bn^{\ell k} \\
&= \det(\gamma) \epsilon_{ij\ell}^{ijm}\gamma^{(d}_m\gamma^{c)}_kn^{\ell k} \\
&= \det(\gamma) \gamma^{(d}_\ell\gamma^{c)}_kn^{\ell k} \\
&= \det(\gamma) \gamma^d_\ell\gamma^c_kn^{\ell k}.\tag{6.1}
\end{align}
Thus, apart from the factor $\det(\gamma)$ we have found that $n^{ij}$ transform as its indices suggest. Note that although the symmetrization vanishes in (6.1) it is precisely the symmetrization that kills the first term in (6). Of course, the transformation of $q_{ij}$ follows trivially from the definition: 
$$
q_{ab} = \gamma^i_a\gamma^j_bq_{ij}.\tag{6.2}
$$
Now note that if we write our quantities as $3\times 3$ matrices we can express (6.1) and (6.2) as
\begin{align}\tag{6.3}
\widetilde{q} &= \gamma q \gamma^T, &
\widetilde{n} &= \det(\gamma) \left(\gamma^{-1}\right)^T n \gamma^{-1},
\end{align}
Under orthogonal transformatios we may rewrite (6.3) as
\begin{align}\tag{6.4}
\widetilde{q} &= \gamma q \gamma^{-1}, &
\widetilde{n} &= \det(\gamma)\gamma n \gamma^{-1},
\end{align}
which, apart from the factor $\det(\gamma)$ is identical to the mutual transformation of two linear operators. Since both $q_{ij}$ and $n^{ij}$ are symmetric they are diagonalized by orthogonal transformations, whence we can confidently say that $q$ and $n$ are mutually diagonalizable if and only if there is a linear transformation that brings them to commute as matrices. Technically we are using the trivial fact that $n^{ab}$ is diagonal if and only if $\det(\gamma)^{-1}n^{ab}$ is.
The trick is then to remember that although $q_{ij}$ and $n^{ij}$ transform as matrices under orthogonal transformations, they do not in general do so. In particular, diagonal transformations allow us to rescale the eigenvalues. More precisely we may assume, without restriction, that $q_{ij}$ has been diagonalized. The commutator components are then given by 
\begin{align}
q_{1i}n^{i2} - n^{1i}q_{i2} &= (q_{11} - q_{22})n^{12}, \\
q_{1i}n^{i3} - n^{1i}q_{i3} &= (q_{11} - q_{33})n^{13}, \\
q_{2i}n^{i3} - n^{2i}q_{i3} &= (q_{22} - q_{33})n^{23},
\end{align}
and under a suitable diagonal transformation (such as a normalization) we may acquire $q_{11} = q_{22} = q_{33}$, thus making the commutator vanish. As noted above, we may then diagonalize $n^{ij}$ without disturbing the orthogonalizaion, and since an orthogonal transformation does not disturb the eigenvalues all the eigenvalues of $q_{ij}$ are identical. We can therefore choose $a = (a,0,0)$ without restriction. We are then free to normalize $n^{ij}$, and if possible $a$ (by (6.1) we can normalize $a$ if and only if either $n^{22} = 0$, $n^{33}=0$, or both), although we may destroy the normalization of our triad in the process (as is also noted by the authors of the article).
