3D. Phase diagram
The phase diagram was obtained by variational means, considering all the states mentioned in Variational analysis of GL functional. Since the 4th order vertex has a pronounced minimum at a particular angle for low enough temperatures, states that can take the full advantage of that dominate. There are only a few states like this: Columnar, Rhombohedral, Cubic, FCC, and Icosahedral (iQC). The first two can adjust their angle to whatever the optimal angle is. The last three are optimal only if $u$ is minimal at the angle $\pi/2$, the icosahedral angle ( approx 63.4 degrees), or the tetrahedral angle (approx 71 degrees).
As discussed above, the stability of the 4th order GL expansion, can be ensured by introducing featureless repulsion constant C.
At $T = 0.05$, introducing local repulsion with $C>10$ makes all the denominators of variational states positive. The icosahedral state is the lowest energy for $Q \approx 1.9 k_F$ up to $C\approx 15$. For larger $C=\lambda_4$, the stripe (smectic) state wins out.
The general feature of the phase diagram that we find is that higher-symmetry states, which in general have larger number of equivalent q components win out for moderate $C$. The examples are both FCC and iQC. They will clearly when their native angles coincide with the angle that minimizes $u(\alpha)$. Further, they have a finite domain of stability around that angle, outside which they either yield to another high-symmetry state, or go to the generic rhombohedral.
All the phase transition are of the first order, expect for the transition between simple cubic and rhombohedral (since cubic is just a special case of rhombohedral). The transition at $T = 0.05$ occurs near $Q = 1.7 k_F$, as can be seen from the plot of the minimizing angle as a function of $Q$.
The simple reasoning for when large number of ordered components is favored is as follows. To minimize GL energy, we need to minimize GL denominator (see Variational analysis of GL functional). For the states with only one inter-q angle $\tilde\alpha$, the denominator is $u(\tilde \alpha) + [u(0) - 2u(\tilde\alpha)]/N$. Clearly, if the second term is positive, it favors large $N$; in the opposite case, $N = 1$. Both iQC and FCC have $N$ which is larger than what is trivially possible in 3D (achieved by rhomboherdal arrangement); thus they are favored at their optimal angles and have local domains of stability around them.
Finally, note that there is no BCC phase in the phase diagram. It would be stabilized by the cubic invariant, had we included it.
Even though BCC phase does not appear explicitly in the phase diagram as the lowest energy phase, it is very close in energy to iQC in the range of $Q$ where the latter is the winner. This may be expected given that they have the same number of components $N$, with four angles $\pi/3$ (nearly the same as iQC), and only one $\pi/2$.
To test the possibility of other states, energy optimization in momentum space was performed, starting from a large discrete set of $\rho_{q_i}$ (for details see Simulations ). In the predicted range of stability of iQC, this method also consistently gives iQC. There is a possibility of distorted iQC states, however, if the minimum on the 4th order vertex does not exactly coincide with the icosahedral angle, see Distortions. That treatment does not include the possible non-coplanar terms that appear e.g. in the FCC case.
( The phase boundaries can be obtained as described in the attached note, if noncoplanar term is not included)