Is a quadratic Majorana Hamiltonian exactly solvable? Given $2N$ Majorana operators $\{a_i\}$ where $i=1,2,3,4,\cdots,2N$
The system Hamiltonian is the most general quadratic form:
$H=\sum A_{ij}a_i a_j$
where $ \{a_i,a_j\}=2\delta_{ij} \quad a^\dagger_i=a_i $ is the Majorana operator comes from a set of fermion operators $c_i \ c^\dagger_i$ and  $ a_i=c^\dagger_i+c_i \quad a_{2i}=i(c^\dagger_i-c_i) $
$A_{ij}$ is chosen to guarantee that $H^\dagger=H$
So, such a Hamiltonian describes a lot of systems.

My question is:
(1)Since it's quadratic, is it exactly solvable? what is people's understanding of it?
By exactly solvable, my intuition is the analogy of the ordinary fermion quadratic problem $H=\sum h_{ij} c^\dagger_i c_j +h.c. $ , after linear transformation we can find a nice basis  $H=\sum \epsilon_n v^\dagger_n v_n +... $  so, $ v^\dagger_n v^\dagger_m \cdots |ground>$ and $<|c^\dagger_i c^\dagger_j \cdots|>$ can be trivially sovled via inverse linear transformation and expansion.
(2)  let's choose $2L$ different Majorana operators, $a_{c_1},a_{c_2},a_{c_3},\cdots,a_{c_{2L}} \in \{a_i\} $ 
I want to calculate ground state average for the product of these operators.
$<G|\prod_{i=1}^{2L} a_{c_i}|G>=\text{complicated function of }\{c_{n};A_{ij}\}$
is this trivial or difficult ? 
Is the "complicated function"  linear function, polynomial, separable, involves determinant, pfaffian ?

I encounter this structure in my quantum Ising model research.
It would be helpful, if someone could show me similar problems as well.
If the problem is  $H=\sum (B_{ij}c^\dagger_i c_j　＋h.c.)$, all the answer seems to be trivial. 
 A: Let's start with your Hamiltonian $i\sum_{jk} A_{jk} a_j a_k$, where I've factored out an $i$ for convenience, and deduce some properties of the matrix $A$. Since the ${a}_{j}$ anticommute, we can assume without loss of generality that $A_{jk}$ is antisymmetric. In fact, if $A$ is any matrix, than the antisymmetric part of $A$, $\frac{1}{2}(A-A^T)$, generates the same Hamiltonian as $A$ up to an irrelevant constant term. We also know $H$ is Hermitian, so that $i\sum_{jk}A_{jk}a_ja_k=(i\sum_{jk}A_{jk}a_ja_k)^\dagger=
-i\sum_{jk}A_{jk}^*a_ka_j=i\sum_{jk}A_{jk}^*a_ja_k$. Thus, we conclude that the matrix $A$ must have all real entries. Therefore, $A$ is a real, antisymmetric matrix.
It is a theorem that real antisymmetric matrices can be put into a near-diagonal form by a real orthogonal transformation. To be precise,
$$
A = O^T \left[\begin{matrix}
0 & \lambda_1 & &\cdots&&&0\\
-\lambda_1 & 0 & & &&&\vdots\\
 & & 0 & \lambda_2 &&&\\
\vdots&&-\lambda_2&0\\
&&&&\ddots\\
&&&&&0&\lambda_n\\
0&\cdots&&&&-\lambda_n&0
\end{matrix}\right]O
$$
Using this decomposition, we can define new Majorana operators $\bar a_i=\sum_j O_{ij}a_j$. Using the fact that $O$ is real orthogonal, you can prove that these $\bar a$ are Majorana operators as well (they square to one, they are their own conjugate, and they anticommute with each other). In terms of these new operators, our Hamiltonian becomes very simple:
$$
H= i\sum_n \lambda_n(\bar a_{2n}\bar a_{2n+1}-\bar a_{2n+1}\bar a_{2n})=2i\sum_n\lambda_n\bar a_{2n}\bar a_{2n+1}
$$
We've thus taken a Hamiltonian that consisted of $2n$ coupled Majorana operators, and transformed it into a Hamiltonian that consists of $n$ decoupled systems, each with two Majorana operators.
We can then transform it back into a Fermion Hamiltonain by writing $\bar{c}_n=\frac{\bar{a}_{2n}+i\bar{a}_{2n+1}}{2}$ and noting that $i\bar a_{2n}\bar a_{2n+1} = 2\bar{c}_n^\dagger\bar{c}_n+\text{constant}$. Thus,
$$
H = \sum_{n}4\lambda_n\bar{c}_n^\dagger\bar{c}_n+\text{constant}
$$
Note that if you want to calculate the ground state expectation value of products of $a$ operators, it's not so bad. You can write your $a$ operators out in terms of the $\bar{a}_i$s, and you can figure out how the $\bar{a}_i$ act on the ground state.
To read more about solving Majorana Hamiltonians, see here: https://arxiv.org/abs/cond-mat/0010440
