Why when applying the Ginzburg-Landau theory for superconductors and expanding the free energy in terms the order parameter $\psi$ one has to consider only the even powers of $|\psi|$?
I suppose it has to do with some symmetry under $\psi\rightarrow -\psi$ of the system (similar to the ferromagnetic case) but I still have not figured out which.
The other commenters have given various reasonable explanations for this. Here is another one:
Since we approximate the free energy as something like a functional version of Taylor expansion in $\psi$, we want the expansion to be differentiable close to the transition (i.e. where the order parameter $\psi$ goes to 0). Note that odd powers of $|\psi|$ are not complex-differentiable functions of $\psi$.
Here's a reference that makes an argument of this nature.
EDIT: Another reference that makes the same argument is Tinkham's book on superconductivity (Chapter 4)
I am not sure that you want more formal explanation but let me sketch. To develop GL expansion by order parameter $\Delta$, one should start from original 4-fermion theory with contact attraction (BCS theory). Then one should perform decoupling in Cooper channel which means Hubbard-Stratonovich transformation, $$\Delta\bar{\psi}_{\uparrow}\bar{\psi}_{\downarrow}+\bar{\Delta}\psi_{\uparrow}\psi_{\downarrow},$$ and $\Delta$ plays role of Cooper pairs field. After decoupling, we can safely integrate out fermion fields by introducing Nambu spinor, $$\Psi=(\psi_{\uparrow}\,\,\bar{\psi}_{\downarrow})^T$$ and its conjugated $\bar{\Psi}$. Integratin over fermion fields gives effective action of theory which is nothing more than GL expansion. Effective action has form $$\mathrm{tr}\ln(1+\mathcal{G}_0^{-1}\Delta),$$ where $$\mathcal{G}_0^{-1}=\begin{pmatrix}-\partial_{\tau}+\partial^2/(2m)+\mu & 0\\ 0 & -\partial_{\tau}-\partial^2/(2m)-\mu\end{pmatrix}.$$ Expanding log term and taking trace, you can easily verify that all odd terms have zero contribution.
To see that effective action coincides with GL expansion, it is worth mentioning that we consider only static configurations, so $\Delta$ does not depend on time and effective action will be $$S_{\text{GL}}=\int d^3x\left[A|\Delta|^2+B|\partial\Delta|^2+C|\Delta|^4\right].$$ Coefficients can be calculated explicitly, gradient term appears if you consider static but spatil inhomogeneous order parameter field.
For me, it is also important that statement about vanishing of odd terms of $\Delta$ in GL expansion does not come from the fact that $n_s\propto |\Delta|^2$. I assume that you derive $n_s$ and I am not sure that quantity $|\Delta|^2$ can be measured.
