Fourth order GL term
Here we derive the form of the 4th order GL term due to the interaction between ions mediated by itinerant electrons.
$$\Delta \Omega^{(4)} = \frac{v_{q_0}^4 \rho_{q_1} \rho_{q_2}\rho_{q_3}\rho_{q_4}}{4 \beta} \sum_{\omega_n, k} G_1(\omega_n) G_2(\omega_n)G_3(\omega_n)G_4(\omega_n) \\\ = \frac{v_{q_0}^4 \rho_{q_1} \rho_{q_2}\rho_{q_3}\rho_{q_4}}{4 \beta} \sum_{\omega_n, k} \frac{1}{(i\omega_n - \epsilon_1 )(i\omega_n - \epsilon_2 )(i\omega_n - \epsilon_3 )(i\omega_n - \epsilon_4 )}\\\ \equiv \mu(\{q_i\})\rho_{q_1} \rho_{q_2}\rho_{q_3}\rho_{q_4},$$
where $\epsilon_1 = \epsilon_k$, $\epsilon_2 = \epsilon_{k - q_1}$, $\epsilon_3 = \epsilon_{k - q_1 - q_2}$, $\epsilon_4 = \epsilon_{k - q_1 - q_2 - q_3}$, and $q_1 + q_2 + q_3 + q_4 = 0$.
Since densities condense only with a preferred wavevector magnitude, $|q_{1,2,3,4} | = q_0$, and sum of $q_i$ is 0, function $\mu$ only depends on $q_0$ and 1 angle in 2D case, and 2 angles in 3D case.
4th order term include self-interaction and mutual interaction between $\rho_q$'s. Self interaction is generated by box diagrams with momentum transfers $(q_1, q_1, -q_1, -q_1)$ (type A) and $(q_1, -q_1, q_1, -q_1)$ (type B).
The combinatorial multiplicities of these diagrams can be calculated as follows. At every vertex of the box diagram we can place $\rho_{\pm q_1}$. For A type, there the sign pattern has to be such that same signs are adjacent, while for B - interlaced. There are 4 ways to place two adjacent ++ in 4 boxes. (++--), (-++-), (--++), (+--+). Hence A type has multiplicity 4. For B type, there are only two distinct ways to arrange: (+-+-), (-+-+), and multiplicity is 2. Therefore, the self interaction goes as
$\sum_i(4A + 2B)|\rho_{q_i}|^4$
Mutual interaction in 2D depends on the relative angle, $2\alpha$, between $q_1$ and $q_2$ (for $\alpha = 0$ we get self-interaction).
There are three distinct contributions to $\Delta \Omega^{(4)}$, which come from the following arrangements or momenta around the box diagram: $\Delta\Omega_1^{(4)}:(q_1, -q_1, q_2, -q_2)$ (type V$_1$), $\Delta\Omega_2^{(4)}:(q_1, -q_1, -q_2, q_2)$ (type V$_2$), and $\Delta\Omega_3^{(4)}: (q_1, q_2, -q_1, -q_2)$ (type D). Now the combinatorial multiplicities.
V$_1$ + V$_2$: 4 ways to place $q_1$, 2 ways to place $-q_1$ next to it (PBC), 2 way to place $\pm q_2$: total 16. Hence there are 8 diagrams of each type.
D: 4 ways to place $q_1$, 2 ways to place $\pm q_2$. Total multiplicity is 8.
Therefore, the mutual interaction term looks like
$\sum_{i\lt j}(8V_1 + 8V_2 + 8D)|\rho_{q_i}|^2|\rho_{q_j}|^2 = \sum_{i\ne j}(4V_1 + 4V_2 + 4D)|\rho_{q_i}|^2|\rho_{q_j}|^2$.
Notice, that in the limit $q_i \to q_j$, $V_1 \to B$, and $(V_2, D)\to A$.
Hence, going back to the original notation in terms of $u(\alpha)$ , we find that the full 4th order GL term is
$$\delta \Omega^{(4)} = \sum_{i\ne j}u(\alpha_{ij})|\rho_{q_i}|^2|\rho_{q_j}|^2 + \frac 12 \sum_{i}u(0)|\rho_{q_i}|^2 $$
where $u(0) = u(\alpha = 0)$. As we show in the attached note, the limit $\alpha \to 0$ is continuous at finite temperature.
In 3D, there is a possibility of a non-coplanar interaction diagrams. They obtain if there are non-trivial quadruplets of $q_1, ..., q_4$ that add up to 0. Such diagrams exist for example for FCC lattice, which has reciprocal BCC. There are 8 BCC reciprocal vectors, which can be split into two distinct quadruplets (tetrahedra). Each has $4! = 24$ multiplicity. Attached Noncoplanar note describes how the counting goes relative to the coplanar quartic terms.
The function $u(\alpha) $ can be obtained numerically. Then it can be used to determine the stability of different crystalline phases variationally.
The frequency summations could be performed with the help of contour integration (The former two diagrams contain double poles, which have to be treated with care),
$$\Delta \Omega^{(4)}_{1,2} = \frac{v_{q_0}^4 |\rho_{q_1}|^2 |\rho_{q_2}|^2}{4 } \sum_{ k} \frac{n_F(\epsilon_1)}{(\epsilon_1 - \epsilon_2) ^2(\epsilon_1 - \epsilon_4)}+ \frac{n_F(\epsilon_4)}{(\epsilon_4 - \epsilon_1) (\epsilon_4 - \epsilon_2)^2} \\ \ \ \ \ \ \ - \frac{n_F(\epsilon_2)}{(\epsilon_2 - \epsilon_1)^2 (\epsilon_2 - \epsilon_4)} - \frac{n_F(\epsilon_2)}{(\epsilon_2 - \epsilon_1)(\epsilon_2 - \epsilon_4)^2} + -\frac{n_F'(\epsilon_2)}{(\epsilon_2 - \epsilon_1)(\epsilon_2 - \epsilon_4)} . $$
and
$$\Delta \Omega^{(4)}_3 = \frac{v_{q_0}^4 |\rho_{q_1}|^2 |\rho_{q_2}|^2}{4 } \sum_{ k} \frac{n_F(\epsilon_1)}{(\epsilon_1 - \epsilon_2) (\epsilon_1 - \epsilon_3)(\epsilon_1 - \epsilon_4)}+ \frac{n_F(\epsilon_2)}{(\epsilon_2 - \epsilon_1) (\epsilon_2 - \epsilon_3)(\epsilon_2 - \epsilon_4)} \\ \ \ \ \ \ \ + \frac{n_F(\epsilon_3)}{(\epsilon_3 - \epsilon_1) (\epsilon_3 - \epsilon_2)(\epsilon_3 - \epsilon_4)} + \frac{n_F(\epsilon_4)}{(\epsilon_4 - \epsilon_1) (\epsilon_4 - \epsilon_2)(\epsilon_4 - \epsilon_3)}. $$
The non-coplanar term is formally similar to $\Delta \Omega^{(4)}_3$. Note that total multiplicity of non-coplanar term, 24, is the same as the total of all coplanar 4th odder terms, 8+8+8 = 24. Interestingly, non-coplanar term can always be chosen attractive, and hence it favors states where it exists. Cubic example is FCC, and given the constraint that all $|q_i| = q_0$, this appears to be the only option. Without this restriction, clearly any distortion along x,y,z, axes of FCC will preserve the the summation of q's into 0.
In the limit of $T\to 0$, these contributions develop singularity at $2k_F \cos \alpha = |q|$, that has a Fano shape, i.e. goes from strongly repulsive to strongly attractive over the window of momenta that corresponds to energy scale of $T$. We will use this strong angular dependence to try to stabilize QC.
The mutual interaction terms are smoothly connected to the self-interaction ones, as can be shown analytically. See attached note with derivation of various terms and their relationships.